testing Hogg's linear algebra beliefs

beliefs.py

import numpy as np
import scipy.linalg as la

# create data matrix and design matrices A (for sinusoids) and B (for nuisances) and M (for both)
y = np.random.normal(size=(100)) # data
A = np.random.normal(size=(2, 100)) # components we care about
B = np.random.normal(size=(10, 100)) # nuisance components
M = np.vstack((B, A))
C = np.diag(0.25 * np.ones(100))
Cinv = np.diag(4.00 * np.ones(100))
print "A", A.shape
print "B", B.shape
print "M", M.shape
print "Cinv", Cinv.shape

# get the amplitudes by the method used by FM15 and Angus et al
MTM = np.dot(M, np.dot(Cinv, M.T))
print "MTM", MTM.shape
MTy = np.dot(M, np.dot(Cinv, y))
print "MTy", MTy.shape

xbar_M = la.solve(MTM, MTy)
print "best-fit components using M matrix"
print xbar_M[-2:]

xvar_M = la.inv(MTM)
print "marginalized covariance using M matrix"
print xvar_M[-2:, -2:]

# now fire up the http://en.wikipedia.org/wiki/Woodbury_matrix_identity
# - It might be unrecognizable because (a) numpy, (b) stupidity,
#   and (c) there is an infinite matrix involved!
# - We are computing [C + B Lambda B.T].inverse, where Lambda is infinity.
BTB = np.dot(B, np.dot(Cinv, B.T))
BTBinvBTCinv = la.solve(BTB, np.dot(B, Cinv))
CplusBBTinv = Cinv - np.dot(np.dot(Cinv, B.T), BTBinvBTCinv)
print "CplusBBTinv", CplusBBTinv.shape
print "Cinv has", np.sum(la.eigvalsh(Cinv) > 1e-9), "non-zero eigenvalues"
print "CplusBBTinv has", np.sum(la.eigvalsh(CplusBBTinv) > 1e-9), "non-zero eigenvalues"

# now get the amplitudes Hogg's way
ATA = np.dot(A, np.dot(CplusBBTinv, A.T))
print "ATA", ATA.shape
ATy = np.dot(A, np.dot(CplusBBTinv, y))
print "ATy", ATy.shape

xbar_A = la.solve(ATA, ATy)
print "best-fit components using A matrix"
print xbar_A

xvar_A = la.inv(ATA)
print "marginalized covariance using A matrix"
print xvar_A

read it (or clone it) and weep, @dfm. You owe me a beer!

and ps. @RuthAngus we can super-speed-up your code!