Accurate Computations with Totally Nonnegative Matrices
Back
Current version posted: March 4, 2013.
What's new: 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.
Function description:
TNEigenvalues
| Syntax: |
z=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.
|
| 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.
|
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. |
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 and this accuracy
will be the subject of an upcoming paper.
|
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=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.
|
| 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, Massachusetts
Institute of Technology
"firstname"@math.mit.edu
|
|