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 JoseJavier Martinez (see below for what it does).
 October 2017: Fixed mexdlasq1.c to properly call LAPACK's DLASQ1 under 64bit 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 JoseJavier Martinez) for the
bidiagonal decomposition of a rectangular BernsteinVandermonde
matrix.
Acknowledgement
This software was developed under the support of NSF Grant DMS0314286 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), 123.
 Plamen Koev
Accurate computations with totally nonnegative matrices
SIAM J. Matrix Anal. Appl. 29 (2007), 731751.
 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 nbyn 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 (i1)st row to the ith.
The procedure is subtractionfree 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 (i1)st and scaling columns i1 and i by y and 1/y,
respectively.
The procedure is subtractionfree 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 BernsteinVandermonde
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
BernsteinVandermonde linear systems" (Linear Algebra and Its
Applications, 2007).

Comments: 
Written and contributed by Ana Marco and JoseJavier Martinez.

TNBDBVR
Syntax: 
A=TNBDBVR(c,m1) 
Description: 
Computes the bidiagonal decomposition of a
RECTANGULAR BernsteinVandermonde
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 JoseJavier Martinez.

TNBDCV
Syntax: 
A=TNBDCV(d,s,c) 
Description: 
Computes the bidiagonal decomposition of a CauchyVandermonde
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 MartinezPena described in the paper
"Factorizations of CauchyVandermonde matrices with one multiple pole,"
Recent research on pure and applied algebra, 8595, Nova Sci. Publ.,
Hauppauge, NY, 2003.

Comments: 
Written and contributed by Ana Marco, JoseJavier Martinez, Juan
Manuel Pena.

TNbidiag
Syntax: 
C=TNbidiag(B) 
Description: 
Perform GolubKahan 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
10by10 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=[x_{i}^{j1+lambda(nj+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 BjorckPereyratype 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 JoseJavier Martinez. 
TNJInverse
Syntax: 
C=TNJInverse(B) 
Description: 
Given B=BD(A), returns the bidiagonal
decomposition of the Jinverse of A: JA^{1}J, where
J=diag((1)^{i}), to HRA. 
Comments: 
JA^{1}J 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 BjorckPereyra 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=PLDL^{T}P^{T} to HRA.

Comments: 
An O(n^{3}) algorithm for computing a rankrevealing
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=[x_{i}^{j1}] 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

