Matrix Computation in the Wolfram Language

Itai Seggev, Wolfram Research

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

ABSTRACT

Mathematica and the Wolfram Language have long had strong capabilities in linear algebra.  In this talk we will survey the scope of these capabilities, including long-time strengths, recent feature additions (including new in 14!) and some future directions we're investigating.

Introduction

Matrix computation and linear algebra have a long history in the Wolfram Language:

Version 1 (1988)

(Pseudo)Inverse, Det/Minors, Transpose

LinearSolve, NullSpace, RowReduce

Eigenvalues, Eigenvectors, Eigensystem

Version 2 (1991)

QRDecomposition

SchurDecomposition

MatrixExp

Version 3 (1996)

JordanDecomposition

LUDecomposition

Version 4 (1999)

Packed arrays!

Tr (?)

Version 5 (2002)

SingularValueDecomposition

More conveniences: CharacteristicPolynomial, MatrixRank, Norm

Sparse arrays!!

Version 6/7/8 (2007—2010)

Matrix generators (e.g., RotationMatrix)

Matrix predicates (e.g, HermitianMatrixQ)

KronckerProduct

MatrixPlot

Version 9 (2012)

MatrixLog, MatrixFunction

StructuredArrays!!!

(SymmetrizedArray)

Two Themes

WL makes linear algebra accessible

Exploiting type / structure can be extremely powerful

Can be done automatically or manually

Outline

More on accessibility

What you have done for me lately?

CenteredInterval

Finite fields

Structured matrices

The undiscovered country

Accessibility

Just being able to perform linear algebra on arbitrary matrix with a single, unified interface was a great advance in 1988.

Over time, we added not only more powerful functions, but also simple functions that made common tasks easier.

In Version 13, a major documentation refresh

Added “insight” examples

Documentation

Created a free, on-line course in linear algebra

https://www.wolfram.com/wolfram-u/courses/mathematics/introduction-to-linear-algebra/

Uses our tools, but standard content

CenteredInterval

CenteredInterval is a representation of a subset of the reals or complexes with a center and bounding “radius”

Output intervals are guaranteed to be a bound on possible answers (not necessarily tight)

CenteredInterval[1,2]-CenteredInterval[1,2]

Exp[CenteredInterval[1+I,2+I]]

Introduced in 13, linear algebra coming in 14

Key functions:

CharacteristicPolynomial

Det

Eigenvalues, Eigenvectors, Eigensystem

Inverse

LinearSolve

LUDecomposition

MatrixExp

MatrixPower

Example

A matrix with centered intervals:

m=(({{, , }, {, , }, {, , }}));

Decompose the matrix:

{λ,v}=Eigensystem[m];

The eigenvalues are centered intervals:

λ

The eigenvectors:

v//MatrixForm

This should be full of zeroes:

MatrixForm[m-Inverse[v].DiagonalMatrix[λ].v]

And it is:

Map[IntervalMemberQ[#,0]&,%,{2}]//MatrixForm

Bonus: bounds on Exp

The following plot shows that CenteredInterval got the the top and right bounds tight for the Exp example, but did not do as well for the left and bottom bounds:

Show[ParametricPlot[ReIm[Exp[({{\[Piecewise], ({{-1+4t, 0<=t<=1}, {3+2I(t-1), 1<=t<=2}, {3+2I-4(t-2), 2<=t<=3}, {-1+2I-2I(t-3), 3<=t<=4}})}})]],{t,0,4}],Graphics[{FaceForm[None],EdgeForm[Thick],Rectangle@@ReIm[["Bounds"]],AbsolutePointSize[7],Point[ReIm[["Center"]]]}],PlotRange->All]

Finite Fields

Introduced in 13.3; linear algebra coming in 14

Key functions:

CharacteristicPolynomial

Det

Inverse

LinearSolve

LUDecomposition

MatrixPower

MatrixRank

NullSpace

RowReduce

Diff vis-a-vis CenteredInterval

+MatrixRank

+NullSpace

+RowReduce

-Eigensystem

-MatrixExp

Example

A matrix of finite field elements:

m=(({{, , }, {, , }, {, , }}));

Get the lower-upper matrix, permutation, and (irrelevant) condition number:

{lu,p,c}=LUDecomposition[m];

Extract the lower triangular part:

MatrixForm[l=LowerTriangularize[lu,-1]+DiagonalMatrix[{1,1,1}]]

Extract the upper triangular part:

MatrixForm[u=UpperTriangularize[lu]]

Confirm the factorization:

l.u===m[[p]]

PS, arguably, we should have done this for the lower part:

LowerTriangularize[lu,-1]+DiagonalMatrix[{,,}]//MatrixForm

Structured Matrices

SymmetrizedArray introduced in V9

QuantityArray introduced in V10.1

Started introducing structured matrices in 13.1

have added more in each release, including 14

New in 14:

HermitianMatrix

SymmetricMatrix

UnitaryMatrix

OrthogonalMatrix

Complete list

Designed to be user-extensible

An example: why do we care?

DiagonalMatrix returns a normal matrix by default for small dimensions:

DiagonalMatrix[{1,2,3}]

For large dimensions, it returns a structured array:

d=DiagonalMatrix[Range[5000]]

This can be controlled using the TargetStructure option:

DiagonalMatrix[{1,2,3},TargetStructure->"Structured"]

n=DiagonalMatrix[Range[5000],TargetStructure->"Dense"];

Inverting the normal matrix takes a while:

Inverse[n];//AbsoluteTiming

For the structured array, it is a “little” faster:

Inverse[d];//AbsoluteTiming

The structured array can report what functions it is optimized for:

d["StructuredAlgorithms"]

GPU Computing

Wolfram Language has long provided a means for GPU computing with CUDALink and OpenCLLink

Later added CUDA neural network training

This worked well, but required much from users

We’re now looking to package this up to be more accessible...

Below, I show some results of an early experiment hooking up LinearSolve to the GPU

Performance comparison, double-precision

Perform operation at double-precision on both CPU and GPU

This was performed on a high-end GPU (RTX 4090) and high-end CPU (2×Xeon 4410T)

On some lower end cards, the GPU never beats the CPU

Row[#,"   "]&@{ListLinePlot[timings,PlotMarkers->{Automatic, 10},PlotLegends->{"CPU", "GPU"},FrameTicks->{{Automatic,Automatic},{{#,1024#}&/@Range[2,numPoints,6],None}},FrameLabel->Map[Style[#,28]&]@{"dimension","time (s)"},Frame->True,LabelStyle->15,ImageSize->500],ListLinePlot[Divide@@timings,PlotMarkers->{Automatic, 10},FrameTicks->{{Automatic,Automatic},{{#,1024#}&/@Range[2,numPoints,6],None}},FrameLabel->Map[Style[#,28]&]@{"dimension","speedup"},Frame->True,PlotRange->All,LabelStyle->15,ImageSize->500]}

Performance comparison, single-precision

CPU computation done at double-precision, as normal

For GPU computation, time includes conversion to single precision prior to sending to GPU and conversion back to double to return the result

Even with more modest cards, get a nice speed bump

Conclusion

Wolfram Language has long history of bringing innovation to matrix computation and linear algebra

Better performance

Greater accessibility

We continue to innovate