# [SciPy-User] cholesky of sparse symmetric banded matrix

josef.pktd@gmai... josef.pktd@gmai...
Thu Aug 23 14:34:41 CDT 2012

```On Thu, Aug 23, 2012 at 8:19 AM,  <josef.pktd@gmail.com> wrote:
> 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

usage, in case anyone is interested:

estimating loc and scale of a distribution with no or fixed shape
parameters using empirical quantiles.

works also in the case when variance is infinite  cauchy, t(2)
and can also use robust estimation if there are outliers
(reference paper is on Weibull)

uses statsmodels

Josef

>
> Josef
>
>>
>> -n
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
```