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 and Generating Functions:
  • 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).
  • The moment generating function of x is a function of the real vector t and is M(t) = E(exp(t^Tx)) = exp(t^Tm + ½t^TSt).  The moments of x are the derivatives of M(t) evaluated at t=0.
  • The cumulant generating function of x is a function of the real vector t and is g(t) = ln E(exp(t^Tx)) = t^Tm + ½t^TSt. The cumulants of x are the derivatives of g(t) evaluated at t=0.
• 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.
• Distribution of Real and Imaginary Parts:
  • E(x^R) = m^R
  • E(x^I) = m^I
  • Cov(x^R) = Cov(x^I) = E((x^R-m^R)(x^R-m^R)^T) =E((x^I-m^I)(x^I-m^I)^T) = ½S^R
• Distribution of Magnitude: We define s=sqrt(diag(S)) to be a positive real-valued vector of standard deviations and •,  ÷ and ^•2 to be elementwise operators
  • E(abs(x)) = ½ pi^½ s• exp(-abs(m÷s)^•2) • ₁F₁(1.5 , 1 ; abs(m÷s)^•2) where ₁F₁(a,b;z)=M(a,b,z) is the Confluent Hypergeometric or Kummer function (hypergeom.m in MATLAB) [R.16]
  • [m=0] E(abs(x)) = ½ pi^½ s
  • [m=0] Cov(abs(x)) = ¼ pi ss^T • ( ₂F₁([-0.5,-0.5] , 1 ; (ABS(S)÷ss^T )^•2) - 1) where ₂F₁(a,b;z) is the Confluent Hypergeometric function (hypergeom.m in MATLAB) [R.16]

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.