[SciPy-User] linear algebra: quadratic forms without linalg.inv

[josef.pktd at gmail.com]

This is just an exercise.

In econometrics (including statsmodels) we have a lot of quadratic
forms that are usually calculate with a matrix inverse. I finally
spend some time to figure out how to do this with cholesky or LU
composition which should be numerically more stable or accurate (and
faster).

"Don't let that INV go past your eyes" (matlab file exchange)

Josef

""" Use cholesky or LU decomposition to calculate quadratic forms

different ways to calculate matrix product B.T * inv(B) * B

Note: calling convention in sparse not consistent, sparse requires
loop over right hand side

Author: josef-pktd
"""

import numpy as np
from scipy import linalg

B = np.ones((3,2)).T
B = np.arange(6).reshape((3,2)).T

print 'using inv'
Ainv = linalg.inv(A)
print np.dot(Ainv, B[:,0])
print np.dot(Ainv, B)
print reduce(np.dot, [B.T, Ainv, B])

print 'using cholesky'
F = linalg.cho_factor(A)
print linalg.cho_solve(F, B[:,0])
print linalg.cho_solve(F, B)
print np.dot(B.T, linalg.cho_solve(F, B))


print 'using lu'
F = linalg.lu_factor(A)
print linalg.lu_solve(F, B[:,0])
print linalg.lu_solve(F, B)
print np.dot(B.T, linalg.lu_solve(F, B))

from scipy import sparse

Asp = sparse.csr_matrix(A)
print 'using sparse symmetric lu'
F = sparse.linalg.splu(A)
print F.solve(B[:,0])
#print F.solve(B) # wrong results but no exception
AiB = np.column_stack([F.solve(Bcol) for Bcol in B.T])
print AiB
print np.dot(B.T, AiB)

#not:
#Bsp = sparse.csr_matrix(B)
#print B.T * F.solve(Bsp)  # argument to solve must be dense array

print 'using sparse lu'
F = sparse.linalg.factorized(A)
print F(B[:,0])
#print F(B) # wrong results but no exception
AiB = np.column_stack([F(Bcol) for Bcol in B.T])
print np.dot(B.T, AiB)