[SciPy-User] cholesky of sparse symmetric banded matrix
Thu Aug 23 07:19:20 CDT 2012
On Wed, Aug 22, 2012 at 10:19 AM, Nathaniel Smith <firstname.lastname@example.org> wrote:
> On Wed, Aug 22, 2012 at 12:16 AM, <email@example.com> wrote:
>> I would like to get a cholesky decomposition of a symmetric banded matrix
>> and multiply it with a dense array
>> I found
>> and this
>> 1) is there a way to do linear algebra (dot multiplication) directly in the
>> "upper diagonal ordered form"?
>> 2) is there an efficient way to go from the "upper diagonal ordered form" to
>> a sparse diagonal matrix (or both ways)?
>> Is there code that uses this and that I can look at for the pattern?
>> my problem is standard linear least squares, where I have an explicit banded
>> form for the (nobs, nobs) weighting matrix
>> X'WX and X'Wy
>> and I need a transformation X2 = W^(0.5) X and y2 = W^(0.5) y
>> so I get X2'X2 and X2'y2
>> (nobs: number of observations, prime is transpose)
>> My first example only has one upper and one lower off-diagonal, so working
>> with dense is wasteful.
> If you only have one off-diagonal, then you may be best off using
> CHOLMOD on the CSC representation:
> Of course this uses GPLed code.
Thanks, which makes it however difficult or impossible to use with statsmodels
I ended up doing the transformation just with numpy: with banded
matrix, I only need to multiply with diagonals and add
something like this for my special case
def chol_vinv_diags(self, diags):
'''diags is list of 2 diagonals
#what's minimum scipy version ?
from scipy.linalg import cholesky_banded
#use that we only have one off-diagonal
band_low = np.concatenate((diags, diags, )).reshape(2,-1)
result = cholesky_banded(band_low, lower=True)
return result, result[:-1]
def whiten(self, x):
special case diags is banded symmetric with 1 off-diagonal
diags = self.vinv_diags()
chdiags = self.chol_vinv_diags(diags)
if x.ndim == 2:
chdiags = [chdiags[:,None], chdiags[:,None]]
#cholesky has only diagonal and one lower off-diagonal
res = x * chdiags
res[:-1] += x[1:] * chdiags
> SciPy-User mailing list
More information about the SciPy-User