[SciPy-User] cholesky of sparse symmetric banded matrix
josef.pktd@gmai...
josef.pktd@gmai...
Thu Aug 23 07:19:20 CDT 2012
On Wed, Aug 22, 2012 at 10:19 AM, Nathaniel Smith <njs@pobox.com> wrote:
> On Wed, Aug 22, 2012 at 12:16 AM, <josef.pktd@gmail.com> wrote:
>> I would like to get a cholesky decomposition of a symmetric banded matrix
>> and multiply it with a dense array
>>
>> I found
>> http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.cholesky_banded.html
>> and this
>> http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.dia_matrix.html
>>
>> either
>> 1) is there a way to do linear algebra (dot multiplication) directly in the
>> "upper diagonal ordered form"?
>> or
>> 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:
> http://packages.python.org/scikits.sparse/cholmod.html
> http://lists.vorpus.org/pipermail/scikits-sparse-discuss/2012-August/000038.html
>
> 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[0], diags[1], [0])).reshape(2,-1)
result = cholesky_banded(band_low, lower=True)
return result[0], result[1][:-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[0][:,None], chdiags[1][:,None]]
#cholesky has only diagonal and one lower off-diagonal
res = x * chdiags[0]
res[:-1] += x[1:] * chdiags[1]
return res
Josef
>
> -n
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
More information about the SciPy-User
mailing list