• We assume without comment below that A, B, C and their sums and inverses are all symmetric matrices.
• First we define C = (A^-1+B^-1)^-1 = A(BA^-1A+BB^-1A)^-1B = A(B+A)^-1B = A(A+B)^-1B = (similarly) B(A+B)^-1A.
• We also define c = C(A^-1a+B^-1b) = B(A+B)^-1A A^-1a+A(A+B)^-1B B^-1b= A(A+B)^-1b + B(A+B)^-1a
• We see that c^TC^-1c = (a^TA^-1+b^TB^-1)CC^-1C(A^-1a+B^-1b)
  = (a^TA^-1+b^TB^-1)C(A^-1a+B^-1b)
  = a^TA^-1CA^-1a + 2a^TA^-1CB^-1b +b^TB^-1CB^-1b
  = a^T(A+B)^-1BA^-1a + 2a^TA^-1A(A+B)^-1BB^-1b +b^T(A+B)^-1AB^-1b substituting for C
  = a^T(A+B)^-1((A+B)A^-1-I) a + 2a^T(A+B)^-1b +b^T(A+B)^-1((A+B)B^-1-I) b
  = a^T(A^-1- (A+B)^-1) a + 2a^T(A+B)^-1b +b^T(B^-1-(A+B)^-1) b
  = a^TA^-1a +b^TB^-1b -(a^T(A+B)^-1a - 2a^T(A+B)^-1b +b^T(A+B)^-1b)
  = a^TA^-1a +b^TB^-1b -(a-b)^T(A+B)^-1(a-b)
• Hence (x-c)^TC^-1(x-c) = x^TC^-1x -2x^TC^-1c + c^TC^-1c
  = x^T(A^-1+B^-1)x -2x^TC^-1C(A^-1a+B^-1b) + a^TA^-1a +b^TB^-1b -(a-b)^T(A+B)^-1(a-b)
  = x^TA^-1x + x^TB^-1x - 2x^TA^-1a - 2x^TB^-1b + a^TA^-1a + b^TB^-1b - (a-b)^T(A+B)^-1(a-b)
  = (x-a)^TA^-1(x-a) + (x-b)^TB^-1(x-b) - (a-b)^T(A+B)^-1(a-b)
• It follows that N(a ; b , A+B) × N(x ; c, C) = det(2 pi (A+B))^-½
• exp((a-b)^T(A+B)^-1(a-b)) × det(2 pi C)^-½ exp((x-c)^TC^-1(x-c))
  = det(4 pi² (A+B)C)^-½ exp((x-c)^TC^-1(x-c) + (a-b)^T(A+B)^-1(a-b)) multiplying determinants and adding exponents
  = det(4 pi² (A+B)A(A+B)^-1B)^-½ exp((x-a)^TA^-1(x-a) + (x-b)^TB^-1(x-b))
  = det(2 pi A)^-½ exp((x-a)^TA^-1(x-a)) × det(2 pi B)^-½ exp((x-b)^TB^-1(x-b))
  = N(x ; a, A) N(x ; b, B)
  

N(x ; a, A)^m = N(0 ; 0, m(2×pi)^m-2A^m-1) × N(x ; a, m^-1A)

We prove this by induction. For m=1, we have N(x ; a, A)¹ = N(0 ; 0, (2×pi)^-1I) × N(x ; a, A) which is true.

• N(x ; a, A)^m = N(x ; a, A) × N(x ; a, A)^m-1
  = N(x ; a, A) × N(0 ; 0, (m-1)(2×pi)^m-3A^m-2) × N(x ; a, (m-1)^-1A)
  = N(0 ; 0, (m-1)(2×pi)^m-3A^m-2) × N(x ; a, A) × N(x ; a, (m-1)^-1A)
  = N(0 ; 0, (m-1)(2×pi)^m-3A^m-2) × N(a ; a, (1+(m-1)^-1)A) × N(x ; c, C) where from [5.1] C = (A^-1+(m-1)A^-1)^-1 = (mA^-1)^-1 = m^-1A and c = C(A^-1a+(m-1)A^-1a) = mCA^-1a = a
• So N(x ; a, A)^m = N(0 ; 0, (m-1)(2×pi)^m-3A^m-2) × N(a ; a, (1+(m-1)^-1)A) × N(x ; a, m^-1A)
  = N(0 ; 0, 2 × pi × (m-1)(2×pi)^m-3A^m-2× (1+(m-1)^-1)A) × N(x ; a, m^-1A)
  = N(0 ; 0, (m-1)(1+(m-1)^-1)(2×pi)^m-2A^m-1) × N(x ; a, m^-1A)
  = N(0 ; 0, m(2×pi)^m-2A^m-1) × N(x ; a, m^-1A)

• From the definition of a complex gaussian we have
  • E(x_i^Rx_k^R) = ½s_i,k^R + m_i^Rm_k^R
  • E(x_i^Ix_k^I) = ½s_i,k^R + m_i^Im_k^I
  • E(x_i^Ix_k^R) = ½s_i,k^I + m_i^Im_k^R
  • E(x_i^Rx_k^I) = -½s_i,k^I + m_i^Rm_k^I
• E(xx^H)_i,k = E(x_ix_k^C) = E((x_i^Rx_k^R + x_i^Ix_k^I) + j(x_i^Ix_k^R - x_i^Rx_k^I)) = s_i,k^R + js_i,k^I + m_im_k^C = (S + mm^H)_i,k
• E(xx^T)_i,k = E(x_ix_k) = E((x_i^Rx_k^R - x_i^Ix_k^I) + j(x_i^Ix_k^R + x_i^Rx_k^I)) = m_im_k = (mm^T)_i,k

• E(x • x^C) = E(diag(xx^H)) = diag(E(xx^H)) = diag(S + mm^H) = diag(S) + m • m^C  [5.3]

• If we define v = diag(S) and w = m • m^C , then E(p)E(p)^T = (v+w)(v+w)^T
• Assume that y = x - m has zero mean and that E(y • y^C) = diag(S) = v
  • E((y • y^C)(y^T • y^H))_i,k = E(((x_i^R)² + (x_i^I)²)((x_k^R)² + (x_k^I)²)) =  E((x_i^Rx_k^R)² + (x_i^Rx_k^I)² + (x_i^Ix_k^R)² + (x_i^Ix_k^I)²)
  • Now from Wick's theorem and [5.3], E((x_i^Rx_k^R)²) = E(x_i^Rx_i^R)E(x_k^Rx_k^R) + 2(E(x_i^Rx_k^R))² = ¼v_iv_k + ½(s_i,k^R)² = (¼vv^T + ½S^R • S^R)_i,k
  • Similarly E((x_i^Ix_k^I)²) = (¼vv^T + ½S^R • S^R)_i,k and E((x_i^Ix_k^R)²) = E((x_i^Ix_k^I)²) = (¼vv^T + ½S^I • S^I)_i,k
  • Hence E((y • y^C)(y^T • y^H))_i,k = 2(¼vv^T + ½S^R • S^R)_i,k + 2(¼vv^T + ½S^I • S^I)_i,k = (vv^T + S • S^C)_i,k
• Since the expected value of any odd power of y_i is zero and, from [5.3], E(yy^T) = 0, we have
  E((x • x^C)(x^T • x^H)) = E(((y+m) • (y+m)^C)((y+m)^T • (y+m)^H))
  = E((y • y^C)(y^T • y^H) + (y • y^C)(m^T • m^H) + (y • m^C)(y^T • m^H) + (y • m^C)(m^T • y^H) + (m • y^C)(y^T • m^H) + (m • y^C)(m^T • y^H) + (m • m^C)(y^T • y^H) + (m • m^C)(m^T • m^H))
  = vv^T + S • S^C + vw^T + 0 + S • (mm^H)^C + S^C • (mm^H) + 0 + wv^T + ww^T
  = S • S^C + 2(mm^H • S^T)^R + (v+w)(v+w)^T

• E(x • x)_k = E(x_k • x_k) = s_k,k + m_k² = (diag(S) + m • m)_k

• Assume that x = y + m with E(y) = 0 and and Cov(y) = S
• Cov(x • x)_i,k = E(x_i²x_k²) - (s_i,i + m_i²)(s_k,k + m_k²) = E((y_i²+2y_im_i+m_i²)(y_k²+2y_km_k+m_k²)) - (s_i,i + m_i²)(s_k,k + m_k²) = s_i,is_k,k + 2s_i,k² + s_i,im_k² + s_k,km_i² + 4s_i,km_im_k + m_i²m_k² - (s_i,i + m_i²)(s_k,k + m_k²) = 2s_i,k (s_i,k+2m_im_k) = (2 S • (S + 2mm^T))_i,k

• p(x | y) = p(x, y) / p(y) = f₁(y) exp( -½[x-p; y-q]^T S^-1 [x-p; y-q] ) where S = [P R; R^T Q]
• From  [3.5] S^-1 = [C^-1, -C^-1RQ^-1; -Q^-1R^TC^-1, Q^-1(I+R^TC^-1RQ^-1)] where C =(P-RQ^-1R^T)
• Hence [x-p; y-q]^T S^-1 [x-p; y-q] = (x-p)^TC^-1(x-p) - 2(x-p)^TC^-1RQ^-1(y-p) + f₂(y) = (x-p-RQ^-1(y-q))^TC^-1(x-p-RQ^-1(y-q)) + f₃(y)
• Hence p(x | y) =  f₄(y) exp( -½(x-p-RQ^-1(y-q))^TC^-1(x-p-RQ^-1(y-q)) ) = f₅(y) N(x, p+RQ^-1(y-q), C) = f₅(y) N(x, p+RQ^-1(y-q), P - RQ^-1R^T) where f₅(y) must equal 1 to give the correct integral.

• Define z = [x; A^Tx] = [I; A^T](x-m). Then z has mean [m; A^Tm] and Cov(z) = E([I; A^T](x-m)(x-m)^T[I, A]) = [I; A^T]S[I, A] = [S SA; A^TS A^TSA]
• From  [5.8] x | A^Tx=b ~ N(x, m+SA(A^TSA)^-1(b-A^Tm), S - SA(A^TSA)^-1A^TS) = N(x, (I-HA^T)m+Hb, (I-HA^T)S)

• E(y) = E(Ax+b) = A E(x) + b = Am+b
• Var(y) = E(yyT) - E(y)E(y)T
  = E(Axx^TA^T + Axb^T + bx^TA^T + bb^T) - (Amm^TA^T + Amb^T + bm^TA^T + bb^T)
  = (A(S+mm^T)A^T + Amb^T + bm^TA^T + bb^T) - (Amm^TA^T + Amb^T + bm^TA^T + bb^T)
  = ASA^T

• |r_i,j|² = s_i,j² (s_i,is_j,j)^-1 <= 1 from  [3.6]

This proof is adapted from [R.1]. Since variances are unaffected by a shift of mean we can assume wlog that p=0 and q=0.  Now set b = a - Q^-1R_i, then for any a we have

• Now note that for a zero-mean scalar z, Var(z) = E(zz^T). Hence
  Var(x_i - a^Ty) = Var(x_i - R_i^TQ^-1y - b^Ty)
  = E((x_i - R_i^TQ^-1y - b^Ty)(x_i - R_i^TQ^-1y - b^Ty)^T)
  = E(x_ix_i^T + R_i^TQ^-1yy^TQ^-1R_i + b^Tyy^Tb - 2x_iy^TQ^-1R_i - 2x_iy^Tb + 2R_i^TQ^-1yy^Tb)
  = diag(P)_i + R_i^TQ^-1QQ^-1R_i + b^TQb - 2R_i^TQ^-1R_i - 2R_i^Tb + 2R_i^TQ^-1Qb
  = diag(P)_i + R_i^TQ^-1R_i + b^TQb - 2R_i^TQ^-1R_i - 2R_i^Tb + 2R_i^Tb
  = diag(P)_i - R_i^TQ^-1R_i + b^TQb

Only the last term depends on b and, since Q is positive semidefinite, this is minimum when b = 0.

This proof is adapted from [R.1]. Since correlations are unaffected by a shift of mean we can assume wlog that p=0 and q=0. For any a we can define c=sqrt(R_i^TQ^-1R_i / a^TQa). From  [5.12]  we know that

• Var(x_i - ca^Ty) >= Var(x_i - R_i^TQ^-1y)
• E(x_ix_i^T + c²a^Tyy^Ta - 2cx_iy^Ta) >= E(x_ix_i^T +  R_i^TQ^-1yy^TQ^-1R_i - 2x_iy^TQ^-1R_i)
• diag(P)_i + c²a^TQa - 2cE(x_iy^Ta) >= diag(P)_i +  R_i^TQ^-1QQ^-1R_i - 2E(x_iy^TQ^-1R_i)
• c²a^TQa - 2cE(x_iy^Ta) >= R_i^TQ^-1R_i - 2E(x_iy^TQ^-1R_i)
• R_i^TQ^-1R_i  - 2cE(x_iy^Ta) >= R_i^TQ^-1R_i - 2E(x_iy^TQ^-1R_i)
• cE(x_iy^Ta) <= E(x_iy^TQ^-1R_i)
• E(x_iy^Ta) / sqrt(a^TQa) <= E(x_iy^TQ^-1R_i) / sqrt(R_i^TQ^-1R_i )
• E(x_iy^Ta) / sqrt( Var(x_i) × a^TQa) <= E(x_iy^TQ^-1R_i) / sqrt( Var(x_i) × R_i^TQ^-1R_i )
• E(x_iy^Ta) / sqrt( Var(x_i) × Var(a^Ty)) <= E(x_iy^TQ^-1R_i) / sqrt( Var(x_i) × Var(R_i^TQ^-1y))

Thus the correlation between x_i and a^Ty is bounded above by the right hand side and attains this bound when a = Q^-1R_i.  The value of the correlation is given by

• E(x_iy^TQ^-1R_i) / sqrt( Var(x_i) × Var(R_i^TQ^-1y))
  = R_i^TQ^-1R_i / sqrt( diag(P)_i × R_i^TQ^-1R_i)
  = sqrt(R_i^TQ^-1R_i /  diag(P)_i)

• 1 - p_ii^-1 det([p_ii , r_i^T; r_i , Q]) det(Q)^-1
• = 1 - p_ii^-1 ((p_ii - r_i^TQ^-1r_i) det(Q)) det(Q)^-1    [3.1]
• = 1 -  (1 - p_ii^-1r_i^TQ^-1r_i) = p_ii^-1r_i^TQ^-1r_i = (FR)_ii  ÷ p_ii = (g_x|y)_i²

• d/dq (ln(p(X)) = d/dm (ln(p(X)) dm/dq + d/dP (ln(p(X)) dP/dS dS/dq where P = S^-1
• d/dm (ln(p(X)) = d/dM (ln(p(X)) dM/dm
  = - ½ d/dM ( tr((X-M)^T S^-1 (X-M))) dM/dm
  = (S^-1(X-M)):^T dM/dm   [Derivative of Trace]
  = (S^-1(X-M)):^T (1_[k#1] ¤ I)    [Derivative of Linear Function]
  = 1_[k#1]^T(X-M)^TS^-1     [Kroneker Product Identity]
• d/dP (ln(p(X)) = ½k  d/dP (ln(det(2pi×P))) - ½  d/dS (tr((X-M)^T P (X-M)))
  = ½kP^-1:^T  - ½  d/dP (tr((X-M)^T P (X-M)))    [Derivative of Determinant]
  = ½kP^-1:^T  - ½  ((X-M)(X-M)^T):^T   [Derivative of Quadratic]
• dP/dS = -(P ¤ P) = -(S^-1 ¤ S^-1)   [Derivative of Inverse]
• d/dP (ln(p(X)) dP/dS = -½( kS:^T  - ½  ((X-M)(X-M)^T):^T)(S^-1 ¤ S^-1)
  = -½(S^-1( kS  - ½  (X-M)(X-M)^T)S^-1):^T
  = -½(( kS^-1  - ½  S^-1(X-M)(X-M)^T)S^-1):^T

• Since v is the sum of k independent identically distributed zero-mean vectors, we can calculate E(vv^T) for k=1 and then multiply the result by k. For k=1, v simplifies to
   v^T = (x-m)^TS^-1 dm/dq + ½(S - (x-m)(x-m)^T ):^T dP/dq where P = S^-1  (see  [5.15])
• When we form the product vv^T, the cross terms are cubic in (x-m) and therefore have zero mean. We deal with the two squared terms in vv^T separately. For the first term
  • k E (((x-m)^TS^-1 dm/dq)^T ((x-m)^TS^-1 dm/dq))
    = k E ((dm/dq)^T S^-1(x-m)(x-m)^TS^-1 dm/dq)
    = k E ((dm/dq)^T S^-1 dm/dq)
• Since the second term has zero mean, we can write
  • ¼k E {(dP/dq)^T (S - (x-m)(x-m)^T ):  (S - (x-m)(x-m)^T ):^T dP/dq}
    = ¼k (dP/dq)^TE { ((x-m)(x-m)^T ):  ((x-m)(x-m)^T):^T }dP/dq - ¼k(dP/dq)^T  S:  S:^T dP/dq
• We define T = TVEC(n,n) so that T S:  = T^T S:  = S^T: =S: since S and T are symmetrical. [see vectorized transpose]
• To calculate the quartic expectation, E {((x-m)(x-m)^T ):  ((x-m)(x-m)^T):^T }, we use Wick's rule take the three possible pairings of the (x-m) terms and, treating x and y as distinct random vectors, we get:
  • E {((x-m)(x-m)^T ):  ((x-m)(x-m)^T):^T }
    = E {((x-m)(x-m)^T ):  ((y-m)(y-m)^T):^T } + E {((x-m)(y-m)^T ):  ((x-m)(y-m)^T):^T } + E {((x-m)(y-m)^T ):  ((y-m)(x-m)^T):^T }
    = S:S:^T + E {((y-m) ¤ (x-m) )((y-m) ¤ (x-m))^T } + E {((y-m) ¤ (x-m) )((y-m) ¤ (x-m))^TT^T } [see kroneker product and vectorized transpose]
    = S:S:^T + E {((y-m)(y-m) ) ¤ ((x-m)(x-m))^T } (I + T^T)
    = S:S:^T + (S ¤ S) (I + T^T)
  • Substituting this into the second term expression gives
    ¼k (dP/dq)^T(S:S:^T + (S ¤ S) (I + T^T))dP/dq - ¼k(dP/dq)^T  S:  S:^T dP/dq
    = ¼k (dP/dq)^T (S ¤ S) (I + T^T)dP/dq   [Since the first and last terms cancel]
    = ½k (dP/dq)^T(S ¤ S) dP/dq    [Since the symmetry of dP means that T^T dP/dq = dP/dq]
• Putting this all together and using dP/dq = dP/dS dS/dq= -(P ¤ P) dS/dq = -(S^-1 ¤ S^-1) dS/dq   [Derivative of Inverse] we obtain
  J = E {vv^T} = k E ((dm/dq)^T S^-1 dm/dq) + ½k (dP/dq)^T(S ¤ S)dP/dq
  = k E ((dm/dq)^T S^-1 dm/dq) + ½k (dS/dq)^T(S^-1 ¤ S^-1)(S ¤ S)(S^-1 ¤ S^-1)dS/dq
  = k E ((dm/dq)^T S^-1 dm/dq) + ½k (dS/dq)^T(S^-1 ¤ S^-1)dS/dq

• g = E (f) means that dg/dq = E(df/dq + f v^T) = E(f v^T) where v^T = d/dq (ln(p(X)).
• For any arbitrary vector a, we have the scalar equation
  a^T dg/dq J^-1 (dg/dq)^T a = E(a^T f v^TJ^-1 (dg/dq)^T a) = Cov(a^T f, a^T dg/dq  J^-1v)
• Hence from the Cauchy-Schwarz inequality
  (a^T dg/dq J^-1 (dg/dq)^T a)² <= Var(a^T f)Var(a^T dg/dq  J^-1v)
  = (a^TCov(f)a) (a^T dg/dq  J^-1Cov(v) J^-1(dg/dq)^T a)
  = (a^TCov(f)a) (a^T dg/dq  J^-1(dg/dq)^T a)   [Since Cov(v) = J ]
• Hence (a^T dg/dq J^-1 (dg/dq)^T a) <= a^TCov(f)a for any a
• So Cov(f) >= dg/dq J^-1 (dg/dq)^T  where >= represents the Loewner partial order.

• Without loss of generality, we assume that m = 0. The characteristic function is f(t) = E(exp(jt^Tx)) = exp(-½t^TSt) where j=sqrt(-1).
• Differentiating E(exp(jt^Tx)) shows that E(prod(x)) = j^-n dⁿf/dt₁dt₂...dt_n evaluated at t=0.
• Expanding the power series f(t) = exp(-½t^TSt) = sum((-½t^TSt)^r/r!; r=0..inf). Term r in this contains only terms in t_x^r where x is some combination of the subscripts 1:n. When differentiated n times, terms with r < n/2 will become zero, whereas terms with r > n/2 will still contain powers of t_x and will therefore be zero when t is set to 0.
• It follows that when n is odd,  all terms vanish and E(prod(x))=0.
• When n is even, E(prod(x)) =  dⁿg/dt₁dt₂...dt_n evaluated at t=0 where g(t) = j^-n (-½t^TSt)^½n/(½n)! = ½^½n(t^TSt)^½n/(½n)!
• We can expand (t^TSt)^½n=t^TSt × t^TSt × ...× t^TSt and we now differentiate in turn with respect to t_n, t_n-1, ... , t₁. We note that t occurs n times in this expansion and that  dt/dt_i = e_i a column of the identity matrix.
• When we differentiate with respect to t_n, we obtain a sum of n terms, in each of which one of the occurences of t has been replaced by e_n and the remaining n-1 occurrences remain as t .
• If we now differentiate with respect to t_n-1, each of the n terms in the previous step changes into a sum of n-1 terms in which one ot the occurences of t is now replaced by e_n-1 and the remaining n-2 occurrences remain as t . We thus have a total of n(n-1) terms.
• We repeat this process for t_n-2, ... , t₁and we end up with n! terms of the form (e_v(1))^TSe_v(2)(e_v(3))^TSe_v(4)...(e_v(n-1))^TSe_v(n) = s_v(1),v(2)s_v(3),v(4)...s_v(n-1),v(n) where v runs through all n! permutations of 1:n.
• Thus E(prod(x)) = (½n)!^-12^-½n sum_v(s_v(1),v(2)s_v(3),v(4)...s_v(n-1),v(n))