Stochastic Matrices

Go to: Introduction, Notation, Index

In all the expressions below, x is a vector of real or complex random variables with whose mean vector and covariance matrix are given by: E(x) = m and Cov(x)=E((x-m)(x-m)^H) = S.

Vectors and matrices a, A, b, B, c, C, d and D are constant (i.e. not dependent on x).

General Properties

• The covariance matrix S is Hermitian and positive semi-definite.
• S is strictly positive definite unless there is a deterministic relation between the elements of x of the form a^Hx = 0 for some non-zero a.
• If the elements of x are uniformly spaced samples from a continuous signal, then S is Toeplitz.
• The symmetric correlation coefficient matrix  (also called correlation matrix) is Corr(x) = DIAG(S)^-½ S DIAG(S)^-½.
  WARNING: Correlation matrix is also used for the matrix E(xx^T) = S + mm^T.
  • The Correlation Coefficient between x_i and x_j equals Corr(x)_i,j = E(x_i x_j)/sqrt(E(x_i²)E(x_i²)) and has magnitude <= 1.  [5.11]
• The precision matrix is T = S^-1.

Special Distributions

The expressions for cubic and quartic expectations given below are restricted to the following special distributions:

Independent

• [x:Real Independent] means that the components of x are real and independent. In particular, we require that E(x_i^px_k^q)=E(x_i^p)E(x_k^q). We define m_r=E((x-m)^r) where the r’th power of the vector is elementwise. Note that S=DIAG(m₂) and m₁=0.

Real Gaussian

• [x:Real Gaussian] means that the components of x_[n] are Real and have a multivariate Gaussian pdf: x ~ N(x ; m, S) =  |2×pi×S|^-½exp( -½(x-m)^T S^-1 (x-m) ) where S is symmetric and +ve semidefinite.
  • ln(p(x)) = -½ ln(det(2pi×S)) - ½(x-m)^T S^-1 (x-m)
• If x is both Gaussian and Independent then m_r = diag( (½S)^½r r! / (½r)!) for even r and 0 for odd r.
• N(x ; m, S) = N(m ; x, S) =  |A| N(Ax+b ; Am+b, ASA^T) for any b and non-singular A.
  • N(x ; m, S) = aⁿ N(ax+b ; am+b, a²S)  for any b and non-zero a.
  • N(x ; m, S) = N(-m ; -x, S) = N(m ; x, S) = N(x+a ; m+a, S) for any a.
• N(Ax; u, R) = N(x; m, S) × N(0; u, R) / N(0; m, S) where S = (A^TR^-1A)^-1 and m = SA^TR^-1u for any A (not necessarily square) with full column rank.
• If x ~ N(x; m, S) then
  • y = Ax+b  ~  N(y; Am+b, ASA^T)   [5.10]
    • y = ax+b  ~  N(y; am+b, a²S)
  • y = F^-T(x-m) ~ N(y; 0, I) where F^TF= S. It is not necessary for F to be symmetric but it can always be chosen to be [see Hermitian].
    • If SQ = QD with Q an orthonormal set of eigenvectors and D diagonal the corresponding positive eigenvalues, then we can define F= D^½Q^T giving F^-T= QD^-½.
  • x | A^Tx=b ~ N(x; (I-HA^T)m+Hb, (I-HA^T)S) where H=SA(A^TSA)⁺ where ()⁺ denotes the pseudoinverse or, if non-singular, the inverse. The symmetry of the covariance may be shown explicitly by writing (I-HA^T)S) = S - SA(A^TSA)⁺A^TS.   [5.9]
    • x | a^Tx=b ~ N(x, m+(b - a^Tm)(a^TSa)^-1×Sa, S - (a^TSa)^-1×Saa^TS)
• Joint Distribution
   If [x; y] ~ N([x; y]; [p; q], [P R^T; R Q]) then in the sections below, we define the regression coefficient matrix of x on y as F=R^TQ^-1 .
  • Linear Sum
    • z = Ax+By+c ~ N(z; Ap+Bq+c, APA^T + BRA^T + AR^TB^T + BQB^T)
  • Conditional Distruibutions
    • x | y ~ N(x; p+F(y-q), P - FR).  [5.8]
      • The mean, p+F(y-q), is the regression function of x on y.
      • The covariance, C = P - FR is is the Schur complement of Q in [P R^T; R Q].
    • y | x ~ N(y; q+RP^-1(x-p), Q - RP^-1R^T). The covariance is the Schur complement of P in [P R^T; R Q].
  • Independence: The following are equivalent:
    1. x and y are independent
    2. R = 0
    3. W = 0 in the precision matrix: T = S^-1 = [P R^T; R Q]^-1 = [U W^T; W V]
    4. N([x; y]; [p; q], [P R^T; R Q]) = N(x; p, P) × N(y; q, Q)
  • Multiple Correlation Coefficients
    The vector of multiple correlation coefficients between x and y has the same dimension as x and is given by g_x|y = sqrt(diag(FR ÷ P) ) where the sqrt() function and  ÷  are elementwise.
    • The minimum (over A) of tr(Cov(x - A^Ty)) is obtained when A = F^T and is equal to tr(P - FR).
    • The maximum (over a) of the correlation between x_i and a^Ty is obtained when a = (F^T)_i = Q^-1r_i and is equal to (g_x|y)_i.   [5.13]
      • All elements of g_x|y lie in the range 0 to 1.
    • diag(Cov(x|y) ÷ Cov(x)) = diag((P - FR) ÷ P) = (1- g_x|y • g_x|y) where  • and  ÷ denote elementwise multiplication and division.
      • Var(x_i|y) = (1-(g_x|y)_i²) Var(x_i) showing that conditioning reduces variance.
    • (g_x|y)_i² = 1 - p_ii^-1 det([p_ii , r_i^T; r_i , Q]) det(Q)^-1  [5.14]
  • Precision Matrix:
    The precision matrix T = S^-1 = [P R^T; R Q]^-1 = [U W^T; W V].
    • If we write the and define, then   [3.5]
      • U = Cov(x | y)^-1 = (P - FR)^-1
      • W = -F^TU
      • V = Q^-1+F^TUF
    • The elements of T are given by
      • t_ii = Var(x_i | x_\i) where x_\i denotes the vector x with x_i deleted.
      • t_ij = -Corr(x_i, x_j | x_\i,j)×(t_ii t_jj)^½
        • t_ij = 0 iff x_i and x_j are conditionally independent given x_\i,j.
    • If I, J and K are a partitioning of the indices 1:n, then
      • The submatrix T_I,J = 0 iff x_I and x_J are conditionally independent: p(x_I, x_J | x_K) = p(x_I | x_K) × p(x_J | x_K)
• Product of Gaussians: N(x ; a, A) N(x ; b, B) = N(a ; b , A+B) × N(x ; c, C) [5.1] where
  • C = (A^-1+B^-1)^-1 = A(A+B)^-1B = B(A+B)^-1A
  • c = C(A^-1a+B^-1b) = A(A+B)^-1b + B(A+B)^-1a
• Power of Gaussian:
  • N(x ; a, A)^m = N(0 ; 0, m(2×pi)^m-2A^m-1) × N(x ; a, m^-1A) [5.2]
    • N(x ; a, A)² = N(0 ; 0, 2A) × N(x ; a, ½A)
• Quotient of Gaussians: N(x ; c, C) / N(x ; a, A) = N(x ; b, B) / N(a ; b , A+B) provided (A-C) is non-singular, where
  • B = (C^-1-A^-1)^-1 = A(A-C)^-1C = C(A-C)^-1A
  • b = B(C^-1c-A^-1a) = A(A-C)^-1c - C(A-C)^-1a
• Characteristic Function: The characteristic function of x is a function of the real vector t and is phi(t) = E(exp(jt^Tx)) = E(cos(t^Tx)) + j E(sin(t^Tx)) = exp(jt^Tm - ½t^TSt) where j=sqrt(-1).
• Differential Entropy
  The differential entropy of x is h(x) =  -E{ln(p(x)} = ½ ln(det(2 pi e S)) nats = ½ log₂(det(2 pi e S)) bits
• Cramer-Rao bound
  Suppose m_[n#1] and S_[n#n] are functions of a parameter vector q_[p#1] and that we take k independent samples of x to form the columns of a data matrix X_[n#k]. In the expressions below,  ¤ denotes the kroneker product, : denotes vectorization and dS/dq is a matrix of dimension n²#p (see derivatives).
  • ln(p(X)) = -½k ln(det(2pi×S)) - ½ tr((X-M)^T S^-1 (X-M))  where  M_[n#k] = m×1_[k#1]^T.
  • The Fisher Score vector, v, is defined by v^T = d/dq (ln(p(X))
    = 1_[k#1]^T(X-M)^TS^-1 dm/dq - ½(k S^-1 - S^-1(X-M)(X-M)^TS^-1):^T dS/dq
    = 1_[k#1]^T(X-M)^TS^-1 dm/dq - ½(k S^-1:^T - ((X-M)(X-M)^T):^T (S^-1 ¤ S^-1)) dS/dq  [5.15]
    • E(v) = 0
    • [k=1] v^T = d/dq (ln(p(X)) =  (x-m)^TS^-1 dm/dq - ½(S^-1 - S^-1(x-m)(x-m)^TS^-1):^T dS/dq
  • The Fisher Information Matrix is defined by J = E(vv^T) = k ((dm/dq)^T S^-1 dm/dq + ½(dS/dq)^T (S ¤ S)^-1 dS/dq ) [5.16]
    • The i,j element of J is given by J_ij =  k ((dm/dq_i)^T S^-1 dm/dq_j + ½ tr(R_iS^-1R_jS^-1)) where R_i satisfies R_i: = dS/dq_i
  • Cramer-Rao bound: If f_[r#1](X) is a function of X with mean value g(q), then Cov(f) >= dg/dq J^-1 (dg/dq)^T where >= represents the Loewner partial order.  [5.17]
    • If g(q) = aq then Cov(f) >= a² J^-1

Complex Gaussian

Definition: In this section, <=> represents the  Complex-to-Real isomporphism and  <-> represents the related vector mapping.

• E(x) = m_[n] <-> a_[2n]
• Cov(x) = E((x-m)(x-m)^H) = E(xx^H) - mm^H = S_[n#n] <=> K_[2n#2n]  [5.3]
• N(y ; a, ½K) = |pi×K|^-½exp( -(y-a)^T K^-1 (y-a) ) = |pi×S|^-1exp( -(x-m)^H S^-1 (x-m) )
  • If S is diagonal (and hence also real) then N(x ; m, S) =  N(y ; a, ½K) =  N(|x-m| ; 0, ½S) × |pi×S|^-½. Thus we can express a complex pdf as a truncated real pdf of the same dimension if the components of x are independent.
• K_[2n#2n] may be divided into 2#2  toeplitz blocks of the form [a -b; b a] (see Givens Rotation)
  • All the 2#2 blocks of K that lie on the main diagonal are positive multiples of I. That is, for each component of x, the real and imaginary parts have the same variance and are uncorrelated.

Linear Expectations

• E(Ax + b) = Am + b
  • E(Ax) = Am
  • E(x + b) = m + b
• Cov(Ax + b) = ASA^T
• E(tr(Y)) = tr(E(Y)) where Y depends on x.

Quadratic Expectations

• E((Ax + a)(Bx + b)^H) = ASB^H + (Am+a)(Bm+b)^H
  • E(xx^H) = S + mm^H
  • E(xa^H x) = (S + mm^H)a
  • E(x^H ax^H) = a^H(S + mm^H)
  • E((Ax)(Ax)^H) = A(S + mm^H)A^H
  • E((x + a)(x + a)^H) = S + (m+a)(m+a)^H
• E((Ax+a)^H (Bx+b)) = tr(ASB^H) + (Am+a)^H (Bm+b)
  • E(x^H x) = tr(S) + m^H m
  • E(x^HAx) = tr(AS) + m^HAm
  • E((Ax)^H (Ax)) = tr(ASA^H) + (Am)^H (Am)
  • E((x+a)^H (x+a)) = tr(S) + (m+a)^H (m+a)
• E((Ax + a) ¤ (Bx + b)) = (A ¤ B) S: + (Am + a) ¤ (Bm + b)
  • E(x ¤ x) = S:

For [x:Real Gaussian] :

• E(x • x) = diag(S) + m • m  = diag(S + m • m^T)    [5.6]
• Cov(x • x) = 2 S • (S + 2mm^T)   [5.7]

For [x:Complex Gaussian] :

• E(xx^T) = mm^T  [5.3]
• E(x • x^C) = diag(S) + m • m^C  [5.4]
• Cov(x • x^C) = E((x • x^C)(x^T • x^H)) - E(x • x^C)E(x • x^C)^T = S • S^C + 2(mm^H • S^T)^R  [5.5]

Cubic Expectations

For [x:Real Independent] :

• E((Ax + a)(Bx + b)^T (Cx + c)) = A DIAG(B^T C) m₃ + tr(BSC^T)×(Am+a) + ASC^T (Bm+b) + (ASB^T +(Am+a)(Bm+b)^T)(Cm+c)
  • E(xx^T x) = m₃ + 2Sm + (tr(S)+ m^T m)×m
  • E((Ax + a)(Ax + a)^T(Ax + a)) = A DIAG(A^T A) m₃ + (2ASA^T + (Am+a)(Am+a)^T)(Am+a) + tr(ASA^T)×(Am+a)
• E((Ax + a)b^T(Cx + c)(Dx + d)^T ) = (Am+a)b^T(CSD^T+(Cm+c) (Dm+d)^T) + (ASC^T+(Am+a)(Cm+c)^T) b (Dm+d)^T + b^T(Cm+c)* (ASD^T - (Am+a)(Dm+d)^T)
  • E(xb^Txx^T) = mb^T(S+mm^T) + (S+mm^T) bm^T + b^Tm* (S - mm^T)

For [x:Real Gaussian] :

• E((Ax + a)(Bx + b)^T(Cx + c)) = ASB^T(Cm+c) + ASC^T(Bm+b) + tr(BSC^T)×(Am+a) + (Am+a)(Bm+b)^T(Cm+c)
  • E(xx^Tx) = 2Sm + (tr(S)+ m^Tm)×m
  • E((Ax + a)(Ax + a)^T(Ax + a)) = (2ASA^T + (Am+a)(Am+a)^T)(Am+a) + tr(ASA^T)×(Am+a)

Quartic Expectations

For [x:Independent] :

For [x:Real Gaussian] :

• E((Ax + a)(Bx + b)^T(Cx + c) (Dx + d)^T) = (ASB^T+(Am+a)(Bm+b)^T)(CSD^T+(Cm+c) (Dm+d)^T) + (ASC^T+(Am+a)(Cm+c)^T)(BSD^T+(Bm+b) (Dm+d)^T) + (Bm+b)^T(Cm+c)×(ASD^T - (Am+a)(Dm+d)^T) + tr(BSC^T)*(ASD^T + (Am+a)(Dm+d)^T)
  • E(xx^Txx^T) = 2(S+mm^T)^2 + m^Tm×(S - mm^T) + tr(S)×(S + mm^T)
  • E(xx^TAxx^T) = E((x^TAx) * xx^T) =(S+mm^T)(A+A^T)(S+mm^T) + m^TAm * (S - mm^T) + tr(AS)×(S + mm^T)
    • E(xx^TAxx^T) = [m=0] SAS + SA^TS + tr(AS)×S
  • E((Ax + a)(Ax + a)^T(Ax + a) (Ax + a)^T) = 2(ASA^T+(Am+a)(Am+a)^T)² + (Am+a)^T(Am+a)×(ASA^T - (Am+a)(Am+a)^T) + tr(ASA^T)×(ASA^T + (Am+a)(Am+a)^T)
• E((Ax + a)^T(Bx + b) (Cx + c)^T(Dx + d)) = tr(AS(C^TD+D^TC)SB^T) + ((Am+a)^TB + (Bm+b)^TA)S(C^T(Dm+d) + D^T(Cm+c)) + (tr(ASB^T)+(Am+a)^T(Bm+b))(tr(CSD^T)+(Cm+c)^T(Dm+d))
  • E(x^Txx^Tx) = 2tr(S²) + 4m^TSm + (tr(S) + m^Tm)²
  • E(x^TAxx^TBx) = tr(AS(B+B^T)S) + m^T(A + A^T)S(B + B^T)m + (tr(AS)+m^TAm)(tr(BS)+m^TBm)
    • E(x^TAxx^TBx) = [m=0]  tr(AS(B+B^T)S) + tr(AS)×tr(BS)
• E(a^Txb^Txc^Txd^Tx) = (a^T(S+mm^T)b)(c^T(S+mm^T)d)+(a^T(S+mm^T)c)(b^T(S+mm^T)d)+(a^T(S+mm^T)d)(b^T(S+mm^T)c)-2a^Tmb^Tmc^Tmd^Tm
• E((Ax + a)^T(Ax + a) (Ax + a)^T(Ax + a)) = 2tr(ASA^TASA^T) + 4(Am+a)^TASA^T(Am+a) + (tr(ASA^T) + (Am+a)^T(Am+a))²

High Powers

For [x:Real Gaussian] :

• [n: odd]  E(prod(x_[n]-m)) = 0.   [5.18]
• [n: even]  E(prod(x_[n]-m)) = (½n)!^-12^-½n sum_v(s_v(1),v(2)s_v(3),v(4)...s_v(n-1),v(n)) where the sum is over all n! permutations v of the numbers 1:n.   [5.18]
  Note that each term in the summation arises (½n)!2^½n times since the ½n factors s_ij can be rearranged in (½n)! orders and for each factor s_ij = s_ji since S is symmetric. Thus an equivalent formula is to omit the normalizing factor,  (½n)!^-12^-½n, and restrict the summation to all distinct pairings of the numbers 1:n. This is Wick's theorem.