4.2 Positive Definite Linear Equations

posv( A, B[, uplo='L'])

Solves

$\begin{displaymath} A X = B, \end{displaymath}$

where A is a real symmetric or complex Hermitian positive definite matrix. On exit, B is replaced by the solution, and A is overwritten with the Cholesky factor. The matrices A and B must have the same type ('d' or 'z'). Raises an ArithmeticError if the matrix is not positive definite.

potrf( A[, uplo='L'])

Cholesky factorization

$\begin{displaymath} A = LL^T \qquad \mbox{or} \qquad A = LL^H \end{displaymath}$

of a positive definite real symmetric or complex Hermitian matrix A. On exit, the lower triangular part of A (if uplo is 'L') or the upper triangular part (if uplo is 'U') is overwritten with the Cholesky factor or its (conjugate) transpose. Raises an ArithmeticError if the matrix is not positive definite.

potrs( A, B[, uplo='L'])

Solves a set of linear equations

$\begin{displaymath} AX=B \end{displaymath}$

with a positive definite real symmetric or complex Hermitian matrix, given the Cholesky factorization computed by posv() or potrf(). On entry, A contains the triangular factor, as computed by posv() or potrf(). On exit, B is replaced by the solution. B must have the same type as A.

potri( A[, uplo='L']): Computes the inverse of a positive definite matrix. On entry, A contains the Cholesky factorization computed by potrf() or posv(). On exit, it contains the inverse.

As an example, we use posv() to solve the linear system

$\begin{displaymath} \left[ \begin{array}{cc} -\mbox{\bf diag}\,(d)^2 & A \\ ... ...ght] = \left[ \begin{array}{c} b_1 \\ b_2 \end{array} \right] \end{displaymath}$

(4.1)

by block-elimination. We first pick a random problem.

>>> from cvxopt.base import matrix, div
>>> from cvxopt.random import normal, uniform
>>> from cvxopt.blas import syrk, gemv
>>> from cvxopt.lapack import posv
>>> m, n = 100, 50  
>>> A = normal(m,n)
>>> b1, b2 = normal(m), normal(n)
>>> d = uniform(m)

We then solve the equations

$\begin{displaymath} A^T \mbox{\bf diag}\,(d)^{-2}A x_2 = b_2 + A^T \mbox{\bf di... ...(d)^{-2} b_1, \qquad \mbox{\bf diag}\,(d)^2 x_1 = Ax_2 - b_1. \end{displaymath}$

>>> Asc = div(A, d[:, n*[0]])                # Asc := diag(d)^{-1}*A
>>> B = matrix(0.0, (n,n))
>>> syrk(Asc, B, trans='T')                  # B := Asc^T * Asc = A^T * diag(d)^{-2} * A
>>> x1 = div(b1, d)                          # x1 := diag(d)^{-1}*b1
>>> x2 = +b2
>>> gemv(Asc, x1, x2, trans='T', beta=1.0)   # x2 := x2 + Asc^T*x1 = b2 + A^T*diag(d)^{-2}*b1 
>>> posv(B, x2)                              # x2 := B^{-1}*x2 = B^{-1}*(b2 + A^T*diag(d)^{-2}*b1)
>>> gemv(Asc, x2, x1, beta=-1.0)             # x1 := Asc*x2 - x1 = diag(d)^{-1} * (A*x2 - b1)
>>> x1 = div(x1, d)                          # x1 := diag(d)^{-1}*x1 = diag(d)^{-2} * (A*x2 - b1)

There are separate routines for equations with positive definite band matrices.

pbsv( A, B[, uplo='L'])

Solves

$\begin{displaymath} AX=B \end{displaymath}$

where A is a real symmetric or complex Hermitian positive definite band matrix. On entry, the diagonals of A are stored in A, using the BLAS format for symmetric or Hermitian band matrices (see section 3.1). On exit, B is replaced by the solution, and A is overwritten with the Cholesky factor (in the BLAS format for triangular band matrices). The matrices A and B must have the same type ('d' or 'z'). Raises an ArithmeticError if the matrix is not positive definite.

pbtrf( A[, uplo='L'])

Cholesky factorization

$\begin{displaymath} A = LL^T \qquad \mbox{or} \qquad A = LL^H \end{displaymath}$

of a positive definite real symmetric or complex Hermitian band matrix A. On entry, the diagonals of A are stored in A, using the BLAS format for symmetric or Hermitian band matrices. On exit, A contains the Cholesky factor, in the BLAS format for triangular band matrices. Raises an ArithmeticError if the matrix is not positive definite.

pbtrs( A, B[, uplo='L'])

Solves a set of linear equations

$\begin{displaymath} AX=B \end{displaymath}$

with a positive definite real symmetric or complex Hermitian band matrix, given the Cholesky factorization computed by pbsv() or pbtrf(). On entry, A contains the triangular factor, as computed by pbsv() or pbtrf(). On exit, B is replaced by the solution. B must have the same type as A.

The following functions are useful for tridiagonal systems.

ptsv( d, e, B)

Solves

$\begin{displaymath} A X = B, \end{displaymath}$

where A is an n by n real symmetric or complex Hermitian tridiagonal matrix, with diagonal d (a 'd' matrix of length n) and subdiagonal e (a 'd' or 'z' matrix of length n-1). The arguments e and B must have the same type. On exit d contains the diagonal elements of D in the LDL ${}\mathrm{^T}$ or LDL ${}\mathrm{^H}$ factorization of A, and e contains the subdiagonal elements of the unit lower bidiagonal matrix L. B is overwritten with the solution X. Raises an ArithmeticError if the matrix is singular.

pttrf( d, e): LDL ${}\mathrm{^T}$ or LDL ${}\mathrm{^H}$ factorization of an n by n real symmetric or complex Hermitian tridiagonal matrix A. On entry, the argument d is a 'd' matrix with the diagonal elements of A. The argument e is 'd' or 'z' matrix with the subdiagonal elements of A. On exit d contains the diagonal elements of D, and e contains the subdiagonal elements of the unit lower bidiagonal matrix L. Raises an ArithmeticError if the matrix is singular.

gttrs( d, e, B[, uplo='L'])

Solves a set of linear equations

$\begin{displaymath} AX=B \end{displaymath}$

where A is an n by n real symmetric or complex Hermitian tridiagonal matrix, given its LDL ${}\mathrm{^T}$ or LDL ${}\mathrm{^H}$ factorization. The argument d is the diagonal of the diagonal matrix D. The argument uplo only matters for complex matrices. If uplo is 'L', then e contains the subdiagonal elements of the unit bidiagonal matrix L. If uplo is 'U', then e contains the complex conjugates of the elements of the unit bidiagonal matrix L. On exit, B is overwritten with the solution X. B must have the same type as e.