Accurate Computations with Totally Nonnegative Matrices of any Rank



Current version posted: June 27, 2022.

Acknowledgement

This software was developed under the support of the San Jose State University Woodward Fund and San Jose State RSCA Funds.

How to cite

  • P. Koev, Accurate eigenvalues and exact zero Jordan blocks of totally nonnegative matrices. Numerische Mathematik, 141: pp 693-713, 2019.

Installation instructions:

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

Description:

This software package can performs matrix computations with potentially singular totally nonnegative matrices to high relative accuracy.
See the above referenced paper from my "Publications" page for a description on how the algorithms work and how to supply the input. All TN matrices are represented by their BIDIAGONAL decompositions. Definitions of the matrices E_i and J_i:
  • The matrix J_i(x,y,z) equals the identity with the only exception of the entries x, y, and z in positions (i,i-1), (i-1,i-1), and (i,i), respectively.
  • E_i(x,y)=J_i(x,y,1);

HRA stands for "High Relative Accuracy." Most functions deliver such accuracy.

Function description:

STNEigenvalues

Syntax: [z,jb]=STNEigenvalues(B,C)
Description: Computes the eigenvalues of the TN matrix A, whose bidiagonal decomposition is stored in [B,C] (i.e., [B,C]=BD(A)) to HRA. z is a vector of eigenvalues, jb is a vector with the sizes of the Jordan blocks corresponding to eigenvalue 0.
Example: [z,jb]=STNEigenvalues([1 0 3; 1 1 0; 1 1 0], [1 1 1; 0 1 1; 1 0 1]) returns z=[3;0;0]; jb=2, meaning that the eigenvalues are 3,0,0 and there is one Jordan block of size 2 for the eigenvalue 0 for the matrix [0 0 0; 1 0 0; 1 1 3], whose bidiagonal decomposition is stored in the matrices passed as arguments.

STNAddToNext

Syntax: [B1,C1]=TNAddToNext(B,C,x,y,i)
Description: If A is a TN matrix, such that [B,C]=BD(A), this function computes the bidiagonal decomposition of E_i(x,y)*A. The procedure is subtraction-free and the result is computed to HRA.
Comments: For the same operation on the columns, work on B',C', and then transpose the result.

STNAddToPrevious

Syntax: [B1,C1]=STNAddToPrevious(B,C,x,y,z,i)
Description: Given [B,C]=BD(A), where A is TN, this routine computes the bidiagonal decomposition of A*J_i(x,y,z). The procedure is subtraction-free and the result is computed to HRA.
Comments: For the same operations on the rows, work on B',C', and then transpose the result.

STNBD

Syntax: [B,C]=STNBD(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. In exact arithmetic, [B,C]=STNBD(A); A1=STNExpand(B,C); will result in A1=A.

STNDiagonalScale

Syntax: [B1,C1]=STNDiagonalScale(B,C,d,i)
Description: If [B,C]=BD(A), computes the bidiagonal decomposition of the matrix resulting from A by scaling column i by a scalar d>=0. Computed to HRA.

STNExpand

Syntax: A=STNExpand(B,C)
Description: Returns the TN matrix A, whose bidiagonal decomposition is stored in [B,C], to HRA.
Comments: Once A is computed, it is not, in general, possible to recover BD(A) from A to HRA.

STNProduct

Syntax: [B,C]=STNProduct(B1,C1,B2,C2)
Description: Given the bidiagonal decompositions [B1,C1]=BD(F) and [B2,C2]=BD(G), computes the bidiagonal decomposition [B,C] of the product FG to HRA.

STNRandom

Syntax: [B,C]=STNRandom(n,k)
Description: Generates a random bidiagonal decomposition of an non TN matrix with k total zeros among the nontrivial entries of the bidiagonal decomposition.

STNRank

Syntax: r=STNRank(B,C)
Description: Computes the rank of the TN matrix whose bidiagonal decomposition is stored in [B,C].

STNSetToZeroAndSwapWithNext

Syntax: [B1,C1]=STNSetToZeroAndSwapWithNext(B,C,i)
Description: Sets the ith column to 0 and swaps it with the following column. It works fine also if the ith column is already 0. Works only for i

STNSetToZeroAndSwapWithPrevious

Syntax: [B1,C1]=STNSetToZeroAndSwapWithPrevious(B,C,i)
Description: Sets the ith row to 0 and swaps it with row i-1. It works fine also if the ith row is already 0.

STNSubmatrix

Syntax: [B1,C1]=STNSubmatrix(B,C,i)
Description: If [B,C]=BD(A), computes the bidiagonal decomposition of the matrix obtained by erasing column i of A and adding a zero column in position n to HRA.
Comments: A zero column in position n is added to maintain the square shape of the matrix. Removing the last column of B1 and the last column of C1 would result in the theoretical bidiagonal decomposition of the resulting n x (n-1) matrix. Once the routines are updated to handle rectangular matrices, the need for a zero last column will be eliminated.

STNTests

Syntax: STNTests
Description: Performs a number of tests on the routines in this package. It is not meant to be exhaustive and no guarantees either way as it relies on comparison with conventional routines and condition numbers come into play, but a useful tool should you decide to modify or improve the routines.

STNtridiag

Syntax: [B1,C1]=STNtridiag(B,C)
Description: Computes the bidiagonal decomposition of a tridiagonal matrix with the same eigenvalues as the original matrix following the procedure of Cryer (1976).

STNBD*

Syntax, e.g.,: [B,C]=STNBDLupas(q,x)
Description: Computes the bidiagonal decomposition of a Lupas matrix. Works similarly for STNBDqBernsteinVandermonde, STNBDhBernsteinVandermonde, and STNBDRationalBernsteinVandermonde.

Plamen Koev
Department of Mathematics, San Jose State University
"firstname" dot "lastname"@sjsu.edu