What’s New in Matrix Algebra Functionality

Daniel Lichtblau, Wolfram Research

#WolframTechConf

ABSTRACT

We will briefly review new matrix and linear algebra functionality that was released in version 14.3. We will then proceed to cover further new functionality and enhancements to existing functions that are scheduled for the next release.

Quick note about functions introduced in 14.3

We introduced several matrix-related functions in 14.3.

CompanionMatrix
JordanMatrix
MatrixPolynomialValue
RangeSpace
EigenvalueDecomposition
PolynomialHermiteDecomposition
PolynomialSmithDecomposition
HermiteReduce
SmithReduce
PolynomialHermiteReduce
MatrixMinimalPolynomial
FrobeniusDecomposition

Several came from the Wolfram Function Repository. Also we updated two or three. Since these are all now in our documentation, and we have limited time here, I will not go into further detail about them.

New in 15.0

Here again, several started out their lives in the WFR.

OrderedSchurDecomposition
JordanReduce
FrobeniusReduce
PfaffianDet
PolarDecomposition
RankDecomposition
LDLDecomposition
BunchKaufmanDecomposition
PopovDecomposition

Also we upgraded a few. These include QRDecomposition, LUDecomposition, MatrixFunction and its close cousins.

I should first note that these are all provisional at this time. I am hopeful that all will be approved though. But please be aware that there might be modest changes.

We will now have a look at several of these new functions.

JordanReduce & FrobeniusReduce

These are smaller-footprint variants on their respective decompositions. Often one does not need the conversion matrix that takes the original matrix to its Frobenius (aka rational) decomposition. Likewise for the similarity matrix that converts to the Jordan normal form. And the normal forms themselves are often less costly in speed and memory to compute.

Jordan form of an integer matrix:

m={{3,7,-22,-6},{7,4,-20,-7},{0,0,-3,0},{6,7,-21,-9}};
JordanReduce[m]

This agrees with the second matrix given by JordanDecomposition:

JordanDecomposition[m][[2]]

Frobenius normal form of an integer matrix:

m=(({{4, -1, -4, 9, -2, -3, 6, 1}, {-5, 1, 0, -2, 2, 2, -3, 0}, {8, -2, -5, 14, -4, -6, 10, 2}, {4, -1, -4, 8, -2, -3, 5, 1}, {7, -1, -4, 15, -4, -6, 10, 1}, {5, -1, -4, 7, -3, -3, 4, 1}, {-3, 1, 4, -7, 1, 2, -5, -1}, {15, -3, -8, 27, -8, -11, 20, 4}}));
f=FrobeniusReduce[m];
MatrixForm[f]

RankDecomposition

The rank decomposition of an m ×n matrix a of rank r is a pair of matrices {s,t} with s having dimensions m × r and t having dimensions r × n, such that a=s.t

This can be quite convenient when the matrix has rank considerably lower than either dimension.

A 3×3 matrix of rank 2:

MatrixRank[m={{3/2,19/4,-16/5},{-17/10,67/10,-93/10},{49/10,-865/100,154/10}}]

So the rank decomposition has a 3×2 and a 2×3:

RankDecomposition[m]

As with many matrix rank-related functions (including the close relative RangeSpace), it takes a Tolerance option.

LDLDecomposition

When a matrix m is Hermitian it can often be decomposed into a form that generalizes the Cholesky decomposition: m=l d (l^T) where l is a lower triangular matrix with ones on the diagonal, d is a diagonal matrix, and (l^T) is the transpose of l. An advantage over CholeskyDecomposition is that now the matrix need not be positive definite.

Create a random symmetric matrix:

mat0=RandomReal[{-1,1},{7,7}];
mat=mat0+Transpose[mat0];

By default the result is given using structured arrays:

{l,d}=LDLDecomposition[mat]

Check the result:

l.d.Transpose[l]==mat

It has negative eigenvalues, hence is not positive definite:

Eigenvalues[mat]

So it has no Cholesky decomposition:

CholeskyDecomposition[mat];

BunchKaufmanDecomposition

The Bunch-Kaufman decomposition also applies to Hermitian matrices. It allows a diagonal with 2×2 blocks, hence a block-diagonal rather than strictly diagonal matrix. Also it computes and uses a permutation matrix p. So the result, {l,d,p}, satisfies p m (p^T)=l d (l^T)

Compute the Bunch-Kaufman decomposition of a real symmetric matrix:

SeedRandom[1234];
n=8;
mat0=RandomReal[{-10,10},{n,n}];
mat=mat0+Transpose[mat0];
{lt,bk,perm}=BunchKaufmanDecomposition[mat]

Check the result:

Max[Abs[lt.bk.Transpose[lt]-perm.mat.Transpose[perm]//Chop]]

PopovDecomposition

This decomposition for a matrix of (univariate) polynomials is analogous to lattice reduction for an integer matrix. Only now “size” is based on polynomial degree rather than integer magnitude.

As this is a decomposition, for a matrix m we get a pair {u,p} where p is the Popov form and u m = p.

An important feature is that the conversion matrix u is invertible over the polynomial ring, that is, it has a constant (as in non-polynomial) determinant. This is referred to as being unimodular.

Compute a Popov decomposition of a 2×3 matrix:

m={{x^3-11,x^4-2x+5,x^4-2x^3+3x-7},{x^2+4x+9,x^3-19,x^3+5x-6}};
{u,p}=PopovDecomposition[m]

Check the matrix equation:

Expand[u.m - p]

Check that u is unimodular:

Det[u]

Upgrades small and large

Small: QRDecomposition now accepts a “Full” specification

In this case it returns a pair {q,r} where q is an orthogonal square matrix regardless of the dimension of the input matrix m. Also r will have the same dimensions as m.

Large: MatrixFunction et al use an improved method for exact numeric inputs

This is based on a Wolfram Function Repository function written by Arturas Acus, using a new method he developed with Adolfas Dargys.

mat={{0,0,0,-2,-2,0,2,-1},{0,-2,4,2,-2,2,5,2},{4,0,0,2,4,-3,-2,6},{-2,-2,2,-2,-1,-4,0,2},{0,-2,2,-1,0,2,-2,-2},{0,-4,1,2,-2,2,-4,0},{0,1,0,0,2,0,-4,0},{-1,0,-2,4,-2,4,-4,-2}};
Timing[mexp1=MatrixExp[mat];]
LeafCount[mexp1]
mexp1//N//Chop

Previously this gave a huge result, and quite slowly (even getting the leaf count or evaluating numerically takes time). This next is run in an older kernel:

Timing[mexp=MatrixExp[mat];]
Timing[LeafCount[mexp]]
Timing[mexp//N//Chop]

Absolutely top secret new functionality

Don’t even whisper about this

As several users have requested, we extended the Arnoldi method for computing select eigenvalues to accept a “black-box” function in lieu of an explicit matrix. It is used for computing matrix-times-vector products. Which means no internal code will attempt even a sparse decomposition on the input matrix. The advantage is that are now able to avoid using considerable memory from such a decomposition when the dimensions are significant.

As the “matrix” is now opaque to us, the dimension is also required in the input.

Quit

dim=30;
SeedRandom[12221];
mat=RandomReal[{-10,10},{dim,dim}];
dotFunF[mat_]:=Function[x,mat.x]
LinearAlgebra`MatrixFreeEigenvalues[{dim,dotFunF[mat]},-2]

Since we actually have this matrix, we will compare to default handling by Eigenvalues:

Eigenvalues[mat,-2]

They are the same up to choice of the second smallest. Which is not bad, for Arnoldi.

I expect this to make its way into the Function Repository once the next version has been released.

Summary

In versions 14.3 and 15.0 we are introducing a bunch of new matrix algebra functions (including a Bunch-something function, no less). We hope some of you find these useful. And there are very real upgrades under the hood to several existing functions.

Acknowledgements

I thank Devendra Kapadia and Roger Germundsson for their extensive and careful work on design specifications for these functions. I thank them and also Itai Seggev for all their assistance in documenting these new functions. And I particularly thank Ying Yin in our QA group. Without her patient and painstaking testing and collation of results, these functions would perhaps still exist. But they would not be much fun to use.

Thank you for your attention

MatrixConditionNumber

This is often requested. It will compute either the exact or approximate condition number based on the 1, 2 or Infinity (default) norm. The exact case is only for exact matrices, and is rarely of interest in practice.

For matrices of approximate numbers, the 2-norm case is done using the singular values and it gives a very accurate approximation. The 1 and Infinity norms are cruder approximations using fairly standard methods from the LAPACK software.

The default:

MatrixConditionNumber[m={{1.6,2.7,3.6},{1.2,3.2,5.2},{3.3,3.4,6.5}}]

1-norm:

MatrixConditionNumber[m,1]

2-norm:

MatrixConditionNumber[m,2]

It agrees with the advertised singular values ratio:

svals=SingularValueList[m];
First[svals]/Last[svals]