cholesky of sparse symmetric banded matrix
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_ba... and this http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.dia_matrix.... 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. Josef
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_ba... and this http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.dia_matrix....
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.... Of course this uses GPLed code. -n
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_ba... and this http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.dia_matrix....
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....
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
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_ba... and this http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.dia_matrix....
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....
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) https://github.com/josef-pkt/statsmodels/commit/f59ecd2f1903b91a58d47f94496c... uses statsmodels Josef
Josef
-n _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
participants (2)
-
josef.pktd＠gmail.com -
Nathaniel Smith