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.
-- 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
