# [SciPy-dev] cholesky decomposition for banded matrices

Tim Leslie tim.leslie at gmail.com
Sun Dec 10 17:56:27 CST 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
>
>
>
```