[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