Accurate Computations with Totally Nonnegative Matrices

Back
Current version posted: August 16, 2013.
What's new:
• TNQR now also returns Q. Contributed by Richard Carini.
• TNBDBVR (contributed by Jose-Javier Martinez) for the bidiagonal decomposition of a rectangular Bernstein-Vandermonde matrix.

Acknowledgement

This software was developed under the support of NSF Grant DMS-0314286

Installation instructions:

Download and unzip the package TNTool.zip. You may need to recompile mexdlasq1.c by typing "mex mexdlasq1.c" at the MATLAB prompt.

Description:

This software package can perform virtually all matrix computations with nonsingular totally nonnegative matrices to high relative accuracy. In other words, despite the typical notorious ill conditioning of these matrices, it is possible to perform virtually all computations with them almost as if no rounding errors occur in the process (meaning that the uncertainty in the output is about the same as that in the input).

Most of the algorithms are described here and here. All TN matrices are represented by their BIDIAGONAL decompositions.
HRA stands for "High Relative Accuracy." Most functions deliver such accuracy.

TNEigenvalues

Syntax: z=TNEigenvalues(B) Computes the eigenvalues of the TN matrix A, whose bidiagonal decomposition is stored in B (i.e., B=BD(A)) to HRA. TNEigenvalues([1 2; 3 4]) returns the eigenvalues: [10.6235; 0.3765] of the TN matrix [1 2; 3 10], whose bidiagonal decomposition is stored as [1 2; 3 4]. TNEigenvalues(ones(n)) returns the eigenvalues of the n-by-n Pascal matrix.

TNSingularValues

Syntax: z=TNSingularValues(B) Computes the singular values of the matrix A, such that B=BD(A), to HRA.

Syntax: A=TNAddToNext(B,x,i) If A is a TN matrix, such that B=BD(A), this function computes the bidiagonal decomposition of the TN matrix obtained from A by adding x times the (i-1)st row to the ith. The procedure is subtraction-free and the result is computed to HRA. For the same operations on the columns, work on B' and transpose.

Syntax: A=TNAddToPrevious(B,x,y,i) Given B=BD(A), where A is TN, this routine computes the bidiagonal decomposition of the TN matrix obtained from A, by adding a multiple x of the ith column to the (i-1)st and scaling columns i-1 and i by y and 1/y, respectively. The procedure is subtraction-free and the result is computed to HRA. For the same operations on the rows, work on B' and transpose.

TNBDBV

Syntax: A=TNBDBV(c) Computes the bidiagonal decomposition of a Bernstein-Vandermonde matrix corresponding to a Lagrange interpolation problem using the Bernstein basis. Here c is the vector containing the interpolation nodes. It is based on the paper "A fast and accurate algorithm for solving Bernstein-Vandermonde linear systems" (Linear Algebra and Its Applications, 2007). Written and contributed by Ana Marco and Jose-Javier Martinez

TNBDBVR

Syntax: A=TNBDBVR(c,m1) Computes the bidiagonal decomposition of a RECTANGULAR Bernstein-Vandermonde matrix. Here c is the vector of the nodes and m1 is the number of columns of the matrix. Written and contributed by Ana Marco and Jose-Javier Martinez

TNBDCV

Syntax: A=TNBDCV(d,s,c) Computes the bidiagonal decomposition of a Cauchy-Vandermonde matrix whose only one pole is multiple. d is the pole. s is the multiplicity of d. c is the vector with the interpolation nodes. It uses the algorithm from Martinez-Pena described in the paper "Factorizations of Cauchy-Vandermonde matrices with one multiple pole," Recent research on pure and applied algebra, 85-95, Nova Sci. Publ., Hauppauge, NY, 2003. Written and contributed by Ana Marco, Jose-Javier Martinez, Juan Manuel Pena.

TNbidiag

Syntax: C=TNbidiag(B) Perform Golub-Kahan bidiagonalization of a TN matrix A, such that B=BD(A). Returns BD(C), where C is bidiagonal, computed to HRA.

TNCauchyBD

Syntax: B=TNCauchyBD(x,y) Computes the bidiagonal decomposition B=BD(C) of the Cauchy matrix C(i,j)=1/(x(i)+y(j)) to HRA. One must have x(1)<...< x(n), y(1)<...< y(n) and x(1)+y(1)>0 in order for C to be TN.

TNCauchyEig

Syntax: z=TNCauchyEig(x,y) Computes the eigenvalues of the Cauchy matrix C(i,j)=1/(x(i)+y(j)) to HRA. One must have x(1)<...< x(n), y(1)<...< y(n) and x(1)+y(1)>0 in order for C to be TN. TNCauchyEig([1:10, 0:9]) returns the eigenvalues of the 10-by-10 Hilbert matrix.

TNCauchySVD

Syntax: s=TNCauchySVD(x,y) Computes the singular values of the Cauchy matrix C(i,j)=1/(x(i)+y(j)) to HRA. One must have x(1)<...< x(n), y(1)<...< y(n) and x(1)+y(1)>0 in order for C to be TN.

TNDiagonalScale

Syntax: A=TNDiagonalScale(f,B) If B=BD(A), computes BD(diag(f)*A) to HRA. Works on rectangular matrices.

TNExpand

Syntax: A=TNExpand(B) Returns the TN matrix A, whose bidiagonal decomposition is stored in B, to HRA. Once A is recovered from B=BD(A), it is not, in general, possible to recover BD(A) from A accurately.

TNGVandBD

Syntax: B=TNGVandBD(x,lambda) Computes the bidiagonal decomposition of the TN generalized Vandermonde matrix G=[xij-1+lambda(n-j+1)] to HRA. In order for G to be TN, we must have 0< x(1)<...< x(n) and real, and lambda(1) >= lambda(2) >= ... >= lambda(n) >= 0 must be integers.

TNGVandEig

Syntax: B=TNGVandEig(x,lambda) Computes the eigenvalues of a TN generalized Vandermonde matrix (see TNGVandBD for details).

TNGVandSVD

Syntax: B=TNGVandSVD(x,lambda) Computes the SVD of a TN generalized Vandermonde matrix (see TNGvandBD for details).

TNGVandSolve

Syntax: y=TNGVandSolve(lambda,x,b) Solves Gy=b, where G is a TN generalized Vandermonde matrix (defined as in TNGVandBD), using a Bjorck-Pereyra-type method. See the paper "The Accurate and Efficient Solution of a Totally Positive Generalized Vandermonde Linear System," here for details.

TNInverse

Syntax: C=TNInverse(B) Given B=BD(A), returns C=BD(A-1) to HRA.

TNJInverse

Syntax: C=TNJInverse(B) Given B=BD(A), returns the bidiagonal decomposition of the J-inverse of A: JA-1J, where J=diag((-1)i), to HRA. JA-1J is TN. TNJInverse(B) is the same as abs(TNInverse(B)).

TNProduct

Syntax: C=TNProduct(A,B) Given A=BD(F) and B=BD(G), computes C=BD(FG) to HRA. When F and G are TN, then so is FG.

TNQR

Syntax: [Q,C]=TNQR(B) Given B=BD(A), with A=QR being the QR decomposition of A with the diagonal of R positive, returns C=BD(R) to HRA, and the matrix Q. The matrix Q is also accurate in the appropriate sense, but each element is not accurate to HRA. The computation of Q added on August 16, 2013 by Richard Carini.

TNQRiter

Syntax: C=TNQRiter(B) Given B=BD(A), returns C=BD(F) to HRA, where A is symmetric and F=Q^TAQ is the result of one step of QR iteration without pivoting. In other words, if A=QR is the QR decomposition of A with the diagonal of R being positive, then F=RQ. The matrix F is TP (type 'help TNQRiter' for details).

TNRandom

Syntax: A=TNRandom(m,n) Generates a random TN matrix as A=TNExpand(rand(m,n)).

TNSchurComplement

Syntax: C=TNSchurComplement(B) Given B=BD(A), computes C=BD(E), where E is obtained from A after one step of Gaussian elimination to HRA. One can then "erase" the first row and/or column using TNSubmatrix in order to obtain the bidiagonal decomposition of the Schur complement only.

TNSolve

Syntax: x=TNSolve(B,b) If A is TN and B=BD(A), solves the linear system Ax=b using backward substitution. Thus for any TN matrix for which BD(A) is available, TNSolve yields a Bjorck-Pereyra type method. Solving a Vandermonde system of equations is done as TNSolve(TNVandBD(x),b), where x are the nodes of the Vandermonde matrix.

TNSubmatrix

Syntax: C=TNSubmatrix(B,i) If A is TN and B=BD(A), computes the bidiagonal decomposition of the submatrix of A obtained by erasing row i of A to HRA. A submatrix of a TN matrix is TN. By applying this function repeatedly, one can obtain the bidiagonal decomposition of any submatrix of a TN matrix to HRA, given the bidiagonal decomposition of the original matrix to start with.

TNTridiagGECP

Syntax: [P,L,D]=TNTridiagGECP(B) If T is symmetric, tridiagonal, and TN, and B=BD(T), this routine computes the LDU decomposition of T resulting from Gaussian elimination with complete pivoting T=PLDLTPT to HRA. An O(n3) algorithm for computing a rank-revealing decomposition of a trigiagonal TN matrix (which is simply a positive s.p.d. matrix).

TNtridiag

Syntax: C=TNtridiag(B) If B=BD(A), and A is TN, computes BD(T) to HRA, where T is a symmetric TN tridiagonal matrix similar to A. Building block of TNEigenvalues.

TNTridiagMinor

Syntax: a=TNTridiagMinor(B,alpha,beta) If B=BD(A), where A is tridiagonal symmetric TN matrix, computes det A(alpha,beta) to HRA. alpha and beta are index sets, which must be of the same length and sorted in increasing order.

TNTridiagSubmatrix

Syntax: C=TNTridiagSubmatrix(B,alpha) If B=BD(A), where A is tridiagonal, symmetric, and TN, computes BD(A(alpha,alpha)), where alpha is an index set.

TNVandBD

Syntax: B=TNVandBD(x) Computes the bidiagonal decomposition of the TN Vandermonde matrix V=[xij-1] to HRA. In order for V to be TN, we must have 0< x(1)<...< x(n) and real.

TNVandEig

Syntax: B=TNVandEig(x) Computes the eigenvalues of a TN Vandermonde matrix (see TNVandBD for details).

TNVandSVD

Syntax: B=TNVandSVD(x) Computes the SVD of a TN Vandermonde matrix (see TNVandBD for details).

TNBD

Syntax: B=TNBD(A) Computes the bidiagonal decomposition of the matrix A by performing Neville elimination. This is the only function in the package that does NOT guarantee high relative accuracy. It is provided as a debugging tool, since one must have TNExpand(TNBD(A))=A in exact arithmetic.

Plamen Koev
Department of Mathematics, Massachusetts Institute of Technology
"firstname"@math.mit.edu