[SciPy-dev] cholesky decomposition for banded matrices
Tim Leslie
tim.leslie at gmail.com
Sun Dec 10 18:56:27 EST 2006
On 11/22/06, Jonathan Taylor <jonathan.taylor at stanford.edu> wrote:
> A week or so ago, I asked about generalized eigenvalue problems for
> banded matrices -- turns out all I needed was a Cholesky decomposition.
>
> I added support for banded cholesky decomposition and solution of banded
> linear systems with Hermitian or symmetric matrices in scipy.linalg with
> some tests. The tests are not as exhaustive as they should be....
>
> Patch is attached.
Hi Jonathan,
This patch worked fine for me, so I've checked it in in r2385.
Cheers,
Tim
>
> -- Jonathan
>
>
> Index: Lib/linalg/info.py
> ===================================================================
> --- Lib/linalg/info.py (revision 2324)
> +++ Lib/linalg/info.py (working copy)
> @@ -7,6 +7,7 @@
> inv --- Find the inverse of a square matrix
> solve --- Solve a linear system of equations
> solve_banded --- Solve a linear system of equations with a banded matrix
> + solveh_banded --- Solve a linear system of equations with a Hermitian or symmetric banded matrix, returning the Cholesky decomposition as well
> det --- Find the determinant of a square matrix
> norm --- matrix and vector norm
> lstsq --- Solve linear least-squares problem
> @@ -27,6 +28,7 @@
> diagsvd --- construct matrix of singular values from output of svd
> orth --- construct orthonormal basis for range of A using svd
> cholesky --- Cholesky decomposition of a matrix
> + cholesky_banded --- Cholesky decomposition of a banded symmetric or Hermitian matrix
> cho_factor --- Cholesky decomposition for use in solving linear system
> cho_solve --- Solve previously factored linear system
> qr --- QR decomposition of a matrix
> Index: Lib/linalg/generic_flapack.pyf
> ===================================================================
> --- Lib/linalg/generic_flapack.pyf (revision 2324)
> +++ Lib/linalg/generic_flapack.pyf (working copy)
> @@ -13,6 +13,113 @@
> python module generic_flapack
> interface
>
> + subroutine <tchar=s,d>pbtrf(lower,n,kd,ab,ldab,info)
> +
> + ! Compute Cholesky decomposition of banded symmetric positive definite
> + ! matrix:
> + ! A = U^T * U, C = U if lower = 0
> + ! A = L * L^T, C = L if lower = 1
> + ! C is triangular matrix of the corresponding Cholesky decomposition.
> +
> + callstatement (*f2py_func)((lower?"L":"U"),&n,&kd,ab,&ldab,&info);
> + callprotoargument char*,int*,int*,<type_in_c>*,int*,int*
> +
> + integer optional,check(shape(ab,0)==ldab),depend(ab) :: ldab=shape(ab,0)
> + integer intent(hide),depend(ab) :: n=shape(ab,1)
> + integer intent(hide),depend(ab) :: kd=shape(ab,0)-1
> + integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
> +
> + <type_in> dimension(ldab,n),intent(in,out,copy,out=c) :: ab
> + integer intent(out) :: info
> +
> + end subroutine <tchar=s,d>pbtrf
> +
> + subroutine <tchar=c,z>pbtrf(lower,n,kd,ab,ldab,info)
> +
> +
> + ! Compute Cholesky decomposition of banded symmetric positive definite
> + ! matrix:
> + ! A = U^H * U, C = U if lower = 0
> + ! A = L * L^H, C = L if lower = 1
> + ! C is triangular matrix of the corresponding Cholesky decomposition.
> +
> + callstatement (*f2py_func)((lower?"L":"U"),&n,&kd,ab,&ldab,&info);
> + callprotoargument char*,int*,int*,<type_in_c>*,int*,int*
> +
> + integer optional,check(shape(ab,0)==ldab),depend(ab) :: ldab=shape(ab,0)
> + integer intent(hide),depend(ab) :: n=shape(ab,1)
> + integer intent(hide),depend(ab) :: kd=shape(ab,0)-1
> + integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
> +
> + <type_in> dimension(ldab,n),intent(in,out,copy,out=c) :: ab
> + integer intent(out) :: info
> +
> + end subroutine <tchar=c,z>pbtrf
> +
> + subroutine <tchar=s,d>pbsv(lower,n,kd,nrhs,ab,ldab,b,ldb,info)
> +
> + !
> + ! Computes the solution to a real system of linear equations
> + ! A * X = B,
> + ! where A is an N-by-N symmetric positive definite band matrix and X
> + ! and B are N-by-NRHS matrices.
> + !
> + ! The Cholesky decomposition is used to factor A as
> + ! A = U**T * U, if lower=1, or
> + ! A = L * L**T, if lower=0
> + ! where U is an upper triangular band matrix, and L is a lower
> + ! triangular band matrix, with the same number of superdiagonals or
> + ! subdiagonals as A. The factored form of A is then used to solve the
> + ! system of equations A * X = B.
> +
> + callstatement (*f2py_func)((lower?"L":"U"),&n,&kd,&nrhs,ab,&ldab,b,&ldb,&info);
> + callprotoargument char*,int*,int*,int*,<type_in_c>*,int*,<type_in_c>*,int*,int*
> +
> + integer optional,check(shape(ab,0)==ldab),depend(ab) :: ldab=shape(ab,0)
> + integer intent(hide),depend(ab) :: n=shape(ab,1)
> + integer intent(hide),depend(ab) :: kd=shape(ab,0)-1
> + integer intent(hide),depend(b) :: ldb=shape(b,1)
> + integer intent(hide),depend(b) :: nrhs=shape(b,0)
> + integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
> +
> + <type_in> dimension(nrhs,ldb),intent(in,out,copy,out=x) :: b
> + <type_in> dimension(ldab,n),intent(in,out,copy,out=c) :: ab
> + integer intent(out) :: info
> +
> + end subroutine <tchar=s,d>pbsv
> +
> + subroutine <tchar=c,z>pbsv(lower,n,kd,nrhs,ab,ldab,b,ldb,info)
> +
> + !
> + ! Computes the solution to a real system of linear equations
> + ! A * X = B,
> + ! where A is an N-by-N Hermitian positive definite band matrix and X
> + ! and B are N-by-NRHS matrices.
> + !
> + ! The Cholesky decomposition is used to factor A as
> + ! A = U**H * U, if lower=1, or
> + ! A = L * L**H, if lower=0
> + ! where U is an upper triangular band matrix, and L is a lower
> + ! triangular band matrix, with the same number of superdiagonals or
> + ! subdiagonals as A. The factored form of A is then used to solve the
> + ! system of equations A * X = B.
> +
> + callstatement (*f2py_func)((lower?"L":"U"),&n,&kd,&nrhs,ab,&ldab,b,&ldb,&info);
> + callprotoargument char*,int*,int*,int*,<type_in_c>*,int*,<type_in_c>*,int*,int*
> +
> + integer optional,check(shape(ab,0)==ldab),depend(ab) :: ldab=shape(ab,0)
> + integer intent(hide),depend(ab) :: n=shape(ab,1)
> + integer intent(hide),depend(ab) :: kd=shape(ab,0)-1
> + integer intent(hide),depend(b) :: ldb=shape(b,1)
> + integer intent(hide),depend(b) :: nrhs=shape(b,0)
> + integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
> +
> + <type_in> dimension(nrhs,ldb),intent(in,out,copy,out=x) :: b
> + <type_in> dimension(ldab,n),intent(in,out,copy,out=c) :: ab
> + integer intent(out) :: info
> +
> + end subroutine <tchar=c,z>pbsv
> +
> subroutine <tchar=s,d,c,z>gebal(scale,permute,n,a,m,lo,hi,pivscale,info)
> !
> ! ba,lo,hi,pivscale,info = gebal(a,scale=0,permute=0,overwrite_a=0)
> Index: Lib/linalg/tests/test_cholesky_banded.py
> ===================================================================
> --- Lib/linalg/tests/test_cholesky_banded.py (revision 0)
> +++ Lib/linalg/tests/test_cholesky_banded.py (revision 0)
> @@ -0,0 +1,208 @@
> +import numpy as N
> +import numpy.random as R
> +import scipy.linalg as L
> +from numpy.testing import *
> +
> +from scipy.linalg import cholesky_banded, solveh_banded
> +
> +def _upper2lower(ub):
> + """
> + Convert upper triangular banded matrix to lower banded form.
> + """
> +
> + lb = N.zeros(ub.shape, ub.dtype)
> + nrow, ncol = ub.shape
> + for i in range(ub.shape[0]):
> + lb[i,0:(ncol-i)] = ub[nrow-1-i,i:ncol]
> + lb[i,(ncol-i):] = ub[nrow-1-i,0:i]
> + return lb
> +
> +def _lower2upper(lb):
> + """
> + Convert upper triangular banded matrix to lower banded form.
> + """
> +
> + ub = N.zeros(lb.shape, lb.dtype)
> + nrow, ncol = lb.shape
> + for i in range(lb.shape[0]):
> + ub[nrow-1-i,i:ncol] = lb[i,0:(ncol-i)]
> + ub[nrow-1-i,0:i] = lb[i,(ncol-i):]
> + return ub
> +
> +def _triangle2unit(tb, lower=0):
> + """
> + Take a banded triangular matrix and return its diagonal and the unit matrix:
> + the banded triangular matrix with 1's on the diagonal.
> + """
> +
> + if lower: d = tb[0]
> + else: d = tb[-1]
> +
> + if lower: return d, tb / d
> + else:
> + l = _upper2lower(tb)
> + return d, _lower2upper(l / d)
> +
> +
> +def band2array(a, lower=0, symmetric=False, hermitian=False):
> + """
> + Take an upper or lower triangular banded matrix and return a matrix using
> + LAPACK storage convention. For testing banded Cholesky decomposition, etc.
> + """
> +
> + n = a.shape[1]
> + r = a.shape[0]
> + _a = 0
> +
> + if not lower:
> + for j in range(r):
> + _b = N.diag(a[r-1-j],k=j)[j:(n+j),j:(n+j)]
> + _a += _b
> + if symmetric and j > 0: _a += _b.T
> + elif hermitian and j > 0: _a += _b.conjugate().T
> + else:
> + for j in range(r):
> + _b = N.diag(a[j],k=j)[0:n,0:n]
> + _a += _b
> + if symmetric and j > 0: _a += _b.T
> + elif hermitian and j > 0: _a += _b.conjugate().T
> + _a = _a.T
> +
> + return _a
> +
> +class test_cholesky(NumpyTestCase):
> +
> + def setUp(self):
> + self.a = N.array([N.linspace(0,0.4,5), [1]*5])
> + self.b = N.array([[1]*5, N.linspace(0.1,0.5,5)])
> + self.b[1,-1] = 0.
> +
> +
> + def test_UPLO(self):
> + u, _ = L.flapack.dpbtrf(self.a)
> + l, _ = L.flapack.dpbtrf(self.b, lower=1)
> + assert_almost_equal(band2array(u), band2array(l, lower=1).T)
> +
> + def test_band2array(self):
> + assert_almost_equal(band2array(self.a, symmetric=True), band2array(self.b, lower=1, symmetric=True))
> +
> + def test_factor1(self):
> + c = band2array(L.flapack.dpbtrf(self.a)[0]).T
> + assert_almost_equal(c, N.linalg.cholesky(band2array(self.a, symmetric=True)))
> +
> + def test_factor2(self):
> + c = band2array(L.flapack.dpbtrf(self.a)[0])
> + a = band2array(self.a, symmetric=True)
> + assert_almost_equal(N.dot(c.T, c), a)
> +
> + def test_factor3(self):
> + c = band2array(L.flapack.dpbtrf(self.b, lower=1)[0], lower=1)
> + b = band2array(self.b, symmetric=True, lower=1)
> + assert_almost_equal(N.dot(c, c.T), b)
> +
> + def test_chol_banded1(self):
> + c = cholesky_banded(self.a)
> + a = band2array(self.a, symmetric=True)
> + _c = band2array(c)
> + assert_almost_equal(N.dot(_c.T, _c), a)
> +
> + def test_chol_banded2(self):
> + c = cholesky_banded(self.b, lower=1)
> + a = band2array(self.a, symmetric=True)
> + _c = band2array(c, lower=1)
> + assert_almost_equal(N.dot(_c, _c.T), a)
> +
> + def test_upper2lower(self):
> + assert_almost_equal(self.b, _upper2lower(self.a))
> + assert_almost_equal(self.b, _upper2lower(_lower2upper(self.b)))
> +
> + def test_lower2upper(self):
> + assert_almost_equal(self.a, _lower2upper(self.b))
> + assert_almost_equal(self.a, _lower2upper(_upper2lower(self.a)))
> +
> +class test_complex_cholesky(test_cholesky):
> +
> + def setUp(self):
> + n = 10
> + self.c = N.array([[20]*5,[1j]*4+[0],[0.01]*3+[0]*2])
> +
> + _c = band2array(self.c, hermitian=True, lower=1)
> + self.b = self.c
> + self.a = _lower2upper(self.b)
> +
> +
> + def test_UPLO(self):
> + u, _ = L.flapack.zpbtrf(self.a)
> + l, _ = L.flapack.zpbtrf(self.b, lower=1)
> + assert_almost_equal(band2array(u), band2array(l, lower=1).T)
> +
> + def test_band2array(self):
> + assert_almost_equal(band2array(self.a, hermitian=True), band2array(self.b, lower=1, hermitian=True).conjugate())
> +
> + def test_factor1(self):
> + c = band2array(L.flapack.zpbtrf(self.a)[0]).T.conjugate()
> + assert_almost_equal(c, N.linalg.cholesky(band2array(self.a, hermitian=True)))
> +
> + def test_factor2(self):
> + c = band2array(L.flapack.zpbtrf(self.a)[0])
> + a = band2array(self.a, hermitian=True)
> + assert_almost_equal(N.dot(c.conjugate().T, c), a)
> +
> + def test_factor3(self):
> + c = band2array(L.flapack.zpbtrf(self.b, lower=1)[0], lower=1)
> + b = band2array(self.b, hermitian=True, lower=1)
> + assert_almost_equal(N.dot(c, c.conjugate().T), b)
> +
> + def test_chol_banded1(self):
> + c = cholesky_banded(self.a)
> + a = band2array(self.a, hermitian=True)
> + _c = band2array(c)
> + assert_almost_equal(N.dot(_c.conjugate().T, _c), a)
> +
> + def test_chol_banded2(self):
> + c = cholesky_banded(self.b, lower=1)
> + b = band2array(self.b, hermitian=True, lower=1)
> + _c = band2array(c, lower=1)
> + assert_almost_equal(N.dot(_c, _c.conjugate().T), b)
> +
> +class test_random_cholesky(test_cholesky):
> + def setUp(self):
> + self.c = R.uniform(low=-1,high=0, size=(3,10))
> + self.c[0] = -N.sum(band2array(self.c, lower=1, symmetric=True), axis=1) + 0.1
> + self.b = self.c
> + self.a = _lower2upper(self.b)
> +
> + def test_unit(self):
> + decomp = cholesky_banded(self.c, lower=1)
> + d, l = _triangle2unit(decomp, lower=1)
> + dd = d**2
> + _l = band2array(l, lower=1)
> + assert_almost_equal(N.dot(_l, N.dot(N.diag(dd), _l.T)),
> + band2array(self.c, lower=1, symmetric=True))
> +
> +
> +class test_sym_solver(NumpyTestCase):
> +
> + def test_solveh(self):
> + c, x = solveh_banded(self.a, N.identity(5))
> + a = band2array(self.a, symmetric=True)
> + assert_almost_equal(x, L.inv(a))
> + _c = band2array(c)
> + assert_almost_equal(N.dot(_c.T, _c), a)
> +
> + def setUp(self):
> + self.a = N.array([N.linspace(0,0.4,5), [1]*5])
> + self.b = N.array([[1]*5, N.linspace(0.1,0.5,5)])
> + self.rhs = N.identity(5)
> +
> + def test_solveU(self):
> + c,x,info = L.flapack.dpbsv(self.a, self.rhs)
> + assert_almost_equal(N.dot(x, band2array(self.a, symmetric=True)), N.identity(5))
> +
> + def test_solveL(self):
> + c,x,info = L.flapack.dpbsv(self.b, self.rhs,lower=1)
> + assert_almost_equal(N.dot(x, band2array(self.a, symmetric=True)), N.identity(5))
> + assert_almost_equal(N.dot(x, band2array(self.b, symmetric=True, lower=1)), N.identity(5))
> +
> +if __name__ == "__main__":
> + NumpyTest().run()
> Index: Lib/linalg/basic.py
> ===================================================================
> --- Lib/linalg/basic.py (revision 2324)
> +++ Lib/linalg/basic.py (working copy)
> @@ -10,7 +10,7 @@
> __all__ = ['solve','inv','det','lstsq','norm','pinv','pinv2',
> 'tri','tril','triu','toeplitz','hankel','lu_solve',
> 'cho_solve','solve_banded','LinAlgError','kron',
> - 'all_mat']
> + 'all_mat', 'cholesky_banded', 'solveh_banded']
>
> #from blas import get_blas_funcs
> from flinalg import get_flinalg_funcs
> @@ -170,7 +170,94 @@
> raise ValueError,\
> 'illegal value in %-th argument of internal gbsv'%(-info)
>
> +def solveh_banded(ab, b, overwrite_ab=0, overwrite_b=0,
> + lower=0):
> + """ solveh_banded(ab, b, overwrite_ab=0, overwrite_b=0) -> c, x
>
> + Solve a linear system of equations a * x = b for x where
> + a is a banded symmetric or Hermitian positive definite
> + matrix stored in lower diagonal ordered form (lower=1)
> +
> + a11 a22 a33 a44 a55 a66
> + a21 a32 a43 a54 a65 *
> + a31 a42 a53 a64 * *
> +
> + or upper diagonal ordered form
> +
> + * * a31 a42 a53 a64
> + * a21 a32 a43 a54 a65
> + a11 a22 a33 a44 a55 a66
> +
> + Inputs:
> +
> + ab -- An N x l
> + b -- An N x nrhs matrix or N vector.
> + overwrite_y - Discard data in y, where y is ab or b.
> + lower - is ab in lower or upper form?
> +
> + Outputs:
> +
> + c: the Cholesky factorization of ab
> + x: the solution to ab * x = b
> +
> + """
> + ab, b = map(asarray_chkfinite,(ab,b))
> +
> + pbsv, = get_lapack_funcs(('pbsv',),(ab,b))
> + c,x,info = pbsv(ab,b,
> + lower=lower,
> + overwrite_ab=overwrite_ab,
> + overwrite_b=overwrite_b)
> + if info==0:
> + return c, x
> + if info>0:
> + raise LinAlgError, "%d-th leading minor not positive definite" % info
> + raise ValueError,\
> + 'illegal value in %d-th argument of internal pbsv'%(-info)
> +
> +def cholesky_banded(ab, overwrite_ab=0, lower=0):
> + """ cholesky_banded(ab, overwrite_ab=0, lower=0) -> c
> +
> + Compute the Cholesky decomposition of a
> + banded symmetric or Hermitian positive definite
> + matrix stored in lower diagonal ordered form (lower=1)
> +
> + a11 a22 a33 a44 a55 a66
> + a21 a32 a43 a54 a65 *
> + a31 a42 a53 a64 * *
> +
> + or upper diagonal ordered form
> +
> + * * a31 a42 a53 a64
> + * a21 a32 a43 a54 a65
> + a11 a22 a33 a44 a55 a66
> +
> + Inputs:
> +
> + ab -- An N x l
> + overwrite_ab - Discard data in ab
> + lower - is ab in lower or upper form?
> +
> + Outputs:
> +
> + c: the Cholesky factorization of ab
> +
> + """
> + ab = asarray_chkfinite(ab)
> +
> + pbtrf, = get_lapack_funcs(('pbtrf',),(ab,))
> + c,info = pbtrf(ab,
> + lower=lower,
> + overwrite_ab=overwrite_ab)
> +
> + if info==0:
> + return c
> + if info>0:
> + raise LinAlgError, "%d-th leading minor not positive definite" % info
> + raise ValueError,\
> + 'illegal value in %d-th argument of internal pbtrf'%(-info)
> +
> +
> # matrix inversion
> def inv(a, overwrite_a=0):
> """ inv(a, overwrite_a=0) -> a_inv
>
>
> _______________________________________________
> Scipy-dev mailing list
> Scipy-dev at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-dev
>
>
>
More information about the SciPy-Dev
mailing list