Accurate Computations with Totally Nonnegative Matrices

Back

Current version posted: Febuary, 2019.
What's new:
  • February 2019:
    • A new function TNGVandBD to compute the bidiagonal decomposition of a Generalized Vandermonde matrix, much faster than before, the old version is still available as TNGVandBD_DK05.
    • Updated TNVandBD to return the decompositions of rectangular Vandermonde matrices (not just square).
  • January 2018: Added the function TNInverseExpand, contributed by Ana Marco and Jose-Javier Martinez (see below for what it does).
  • October 2017: Fixed mexdlasq1.c to properly call LAPACK's DLASQ1 under 64-bit architectures. See the notes in the file to properly compile on your system.
  • September 2015: TNtridiag and TNEigenValues now return the transformation and eigenvector matrices, respectively. Research funded by SJSU's Woodward Fund.
  • 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 and the San Jose State University Woodward Fund.

How to cite

  • Plamen Koev, Accurate eigenvalues and SVDs of totally nonnegative matrices SIAM J. Matrix Anal. Appl. 27 (2005), 1-23.
  • Plamen Koev Accurate computations with totally nonnegative matrices SIAM J. Matrix Anal. Appl. 29 (2007), 731-751.
  • The publication referred to in the corresponding routine.

Installation instructions:

Download and unzip the package TNTool.zip. You may need to recompile mexdlasq1.c by typing "mex mexdlasq1.c -lmwlapack" 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.

Function description:

TNEigenvalues

Syntax: [z,V]=TNEigenvalues(B)
Description: Computes the eigenvalues of the TN matrix A, whose bidiagonal decomposition is stored in B (i.e., B=BD(A)) to HRA. It also returns the eigenvectors matrix V, but no claims are made to its accuracy, even though it is computed as a product of highly accurate factors.
Example 1: 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].
Example 2: TNEigenvalues(ones(n)) returns the eigenvalues of the n-by-n Pascal matrix.

TNSingularValues

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

TNAddToNext

Syntax: A=TNAddToNext(B,x,i)
Description: 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.
Comments: For the same operations on the columns, work on B' and transpose.

TNAddToPrevious

Syntax: A=TNAddToPrevious(B,x,y,i)
Description: 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.
Comments: For the same operations on the rows, work on B' and transpose.

TNBDBV

Syntax: A=TNBDBV(c)
Description: 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).
Comments: Written and contributed by Ana Marco and Jose-Javier Martinez.

TNBDBVR

Syntax: A=TNBDBVR(c,m1)
Description: 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.
Comments: Written and contributed by Ana Marco and Jose-Javier Martinez.

TNBDCV

Syntax: A=TNBDCV(d,s,c)
Description: 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.
Comments: Written and contributed by Ana Marco, Jose-Javier Martinez, Juan Manuel Pena.

TNbidiag

Syntax: C=TNbidiag(B)
Description: 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)
Description: Computes the bidiagonal decomposition B=BD(C) of the Cauchy matrix C(i,j)=1/(x(i)+y(j)) to HRA.
Comments: 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)
Description: Computes the eigenvalues of the Cauchy matrix C(i,j)=1/(x(i)+y(j)) to HRA.
Comments: One must have x(1)<...< x(n), y(1)<...< y(n) and x(1)+y(1)>0 in order for C to be TN.
Example: TNCauchyEig([1:10, 0:9]) returns the eigenvalues of the 10-by-10 Hilbert matrix.

TNCauchySVD

Syntax: s=TNCauchySVD(x,y)
Description: Computes the singular values of the Cauchy matrix C(i,j)=1/(x(i)+y(j)) to HRA.
Comments: 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)
Description: If B=BD(A), computes BD(diag(f)*A) to HRA.
Comments: Works on rectangular matrices.

TNExpand

Syntax: A=TNExpand(B)
Description: Returns the TN matrix A, whose bidiagonal decomposition is stored in B, to HRA.
Comments: 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)
Description: 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. Implements the method of Plamen Koev, Accurate computations with totally nonnegative matrices, SIMAX, 2007.

TNGVandBD_DK05

Syntax: B=TNGVandBD(x,lambda)
Description: Same as TNGVandBD above, except it implements an older and slower (but just as accurate) method from Demmel and Koev, The accurate solution to a totally positive generalized Vandermonde linear system, SIMAX, 2005. Provided merely as an implementation of the method in this 2005 paper. The above version, TNGVandBD is better in every way we know of.

TNGVandEig

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

TNGVandSVD

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

TNGVandSolve

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

TNInverse

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

TNInverseExpand

Syntax: C=TNInverseExpand(B)
Description: Given B=BD(A), returns A-1 to HRA. While mathematically equivalent to C=TNExpand(TNInverse(B)), TNInverseExpand is much faster since it produces the inverse directly without having to compute its bidiagonal decomposition first.
Comment: Written and contributed by Ana Marco and Jose-Javier Martinez.

TNJInverse

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

TNProduct

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

TNQR

Syntax: [Q,C]=TNQR(B)
Description: 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.
Notes: The computation of Q added on August 16, 2013 by Richard Carini.

TNQRiter

Syntax: C=TNQRiter(B)
Description: 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)
Description: Generates a random TN matrix as A=TNExpand(rand(m,n)).

TNSchurComplement

Syntax: C=TNSchurComplement(B)
Description: Given B=BD(A), computes C=BD(E), where E is obtained from A after one step of Gaussian elimination to HRA.
Comments: 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)
Description: 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.
Example: 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)
Description: 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.
Comments: 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)
Description: 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.
Comments: 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,V]=TNtridiag(B)
Description: If B=BD(A), and A is TN, computes BD(T) to HRA, where T is a symmetric TN tridiagonal matrix similar to A. It also outputs V, which is such that V^{1}AV=T.
Comments: Building block of TNEigenvalues.

TNTridiagMinor

Syntax: a=TNTridiagMinor(B,alpha,beta)
Description: 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)
Description: 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)
Description: 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)
Description: Computes the eigenvalues of a TN Vandermonde matrix (see TNVandBD for details).

TNVandSVD

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

TNBD

Syntax: B=TNBD(A)
Description: Computes the bidiagonal decomposition of the matrix A by performing Neville elimination.
Comments: 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, San Jose State University
"firstname" dot "lastname"@sjsu.edu