Structured Arrays

Jan Mangaldan

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

ABSTRACT

We discuss recent features in the Wolfram Language related to the representation of and computation with structured arrays. Representing matrices as structured arrays, whenever feasible, often leads to efficiencies in storage, computation time or both. We also touch on their applications.

Introduction

– Structured arrays are a recently added feature in the Wolfram Language (version 13.1+)

Compare with a conventional list:

Introduction

– Preceded by SparseArray (version 5) and SymmetrizedArray (version 9)

IdentityMatrix[5,TargetStructure->SymmetrizedArray]

Introduction

– In general, structured arrays can be generated through an appropriate setting of the TargetStructure option (new in 13.3)

?TargetStructure

– List all the matrices that currently support the TargetStructure option:

Introduction

– Structured arrays are treated the same way as lists when given to other functions

fm=FourierMatrix[128,TargetStructure->"Structured"]

MatrixPlot[fm]

Why have structured arrays?

– They generally take less storage than ordinary lists, but are equivalent to them

– Linear algebra operations on structured arrays generally take less time, memory and effort, due to specialized algorithms implemented within.

Why have structured arrays?

– A structured array:

imStructured=IdentityMatrix[1000,TargetStructure->"Structured"]

– Normal converts structured arrays to lists:

Short[imList=Normal[imStructured]]

– Compare the amount of storage required using ByteCount:

{imStructured,imList}//Map[UnitConvert[N[ByteCount[#]]((Quantity[1, "B", byte, "Bytes"])),"Kilobytes"]&]

Why have structured arrays?

– Compare the time taken for computing the determinant for each:

AbsoluteTiming[Det[imStructured]]

AbsoluteTiming[Det[imList]]

– Compare the time taken for computing the inverse for each:

AbsoluteTiming[imStructuredInv=Inverse[imStructured];]

AbsoluteTiming[imListInv=Inverse[imList];]

Why have structured arrays?

– Structured arrays have properties associated with them:

pm=PermutationMatrix[Cycles[{{1,6,5},{2,3,4}}]]

pm["Properties"]

– Some properties are unique to a particular structured array:

{pm["PermutationCycles"],pm["PermutationList"]}

Why have structured arrays?

– List all built-in functions that will treat the structured array with specialized algorithms using the “StructuredAlgorithms” property:

pm["StructuredAlgorithms"]

– Algorithms on structured arrays will try to return a structured result, whenever applicable:

KroneckerProduct[pm,Transpose[pm]]

Examples

Vandermonde matrix

– Vandermonde matrices take the following form:

VandermondeMatrix[Table[(Subscript[x, k]),{k,5}]]//MatrixForm

– Vandermonde matrices often appear in applications related to interpolating polynomials:

data={{0,0},{0.1,.3},{0.5,.6},{1,-.2},{2,3}};
coefs=LinearSolve[VandermondeMatrix[data[[All,1]]],data[[All,2]]];
poly[x_]=coefs.(x^(Range[0,Length[data]-1]));
Show[Plot[poly[x],{x,0,2},],
ListPlot[data,],]

poly[x]-InterpolatingPolynomial[data,x]//Simplify//Chop

Examples

Vandermonde matrix

– Vandermonde matrices only need the nodes for their construction:

nodes=RandomReal[{-1,1},100];
vmStructured=VandermondeMatrix[nodes]

– Represented as a structured array, the Vandermonde matrix needs little storage compared to the representation as a list:

UnitConvert[Quantity[ByteCount[vmStructured]//N,"Bytes"],"Kilobytes"]

Short[vmList=Normal[vmStructured]]

UnitConvert[Quantity[ByteCount[vmList]//N,"Bytes"],"Kilobytes"]

Examples

Vandermonde matrix

– Operating on Vandermonde matrices as structured arrays usually take less computation time and effort:

AbsoluteTiming[Transpose[vmStructured];]

AbsoluteTiming[Transpose[vmList];]

vec=RandomReal[{-1,1},100];

AbsoluteTiming[LinearSolve[vmStructured,vec];]

AbsoluteTiming[LinearSolve[vmList,vec];]

– The solution of a (dual) Vandermonde system takes O((n^2)) effort (Björck and Pereyra), compared to the usual O((n^3)) effort for solving general linear systems (Gaussian elimination).

Examples

Fourier matrix

– Fourier matrices are matrices with complex exponential elements that appear in the computation of the discrete Fourier transform (DFT):

fm=FourierMatrix[8,TargetStructure->"Structured"]

– Multiplying structured Fourier matrices with a vector is fast because the fast Fourier transform (FFT) is used internally (O(n log(n)) effort):

fmLarge=FourierMatrix[1024,TargetStructure->"Structured"];vec=RandomComplex[{-1-I,1+I},1024];AbsoluteTiming[fmLarge.vec;]

– Matrix-vector multiplication takes O((n^2)) effort in the general case, so the exploitation of structure contributes to efficiency.

Examples

Fourier matrix

– The idea behind FFT is that a Fourier matrix whose dimensions are composite can be factored in terms of Fourier matrices with prime dimensions:

fm2=FourierMatrix[2,TargetStructure->"Structured"]

(factors={KroneckerProduct[fm2,IdentityMatrix[2]],DiagonalMatrix[Flatten[Table[Exp[((2 π I j k)/(2 2))],{k,0,1},{j,0,1}]]],BlockDiagonalMatrix[ConstantArray[fm2,2]],PermutationMatrix[Flatten[Partition[Range[2 2],2],{{2,1}}]]})//Map[MatrixForm]

Dot@@factors==FourierMatrix[2 2]//Simplify

– FFT is effectively equivalent to successively multiplying a vector to be transformed with these structured factors (Van Loan).

Examples

Orthogonal matrix

– Orthogonal matrices can be generated in multiple ways.

– A direct construction from a list of lists:

OrthogonalMatrix[(({{(1/2), (1/2), (1/2), (1/2)}, {(1/2), (1/2), -(1/2), -(1/2)}, {(1/2), -(1/2), (1/2), -(1/2)}, {(1/2), -(1/2), -(1/2), (1/2)}}))]

– From other built-in matrix functions (e.g. rotation and reflection matrices):

RotationMatrix[π/4,{1,-1,1},TargetStructure->"Orthogonal"]

ReflectionMatrix[{1,0,1},TargetStructure->"Orthogonal"]

EulerMatrix[{((3π)/4),(π/3),(π/3)},TargetStructure->"Orthogonal"]

– From matrix decompositions (QR, SVD, etc.):

mat=RandomReal[{-1,1},{4,4}];

QRDecomposition[mat,Pivoting->True,TargetStructure->"Structured"]

SingularValueDecomposition[mat,TargetStructure->"Structured"]

Examples

Orthogonal matrix

– Random orthogonal matrices sampled from the Haar distribution can be generated using CircularRealMatrixDistribution:

om=OrthogonalMatrix[RandomVariate[CircularRealMatrixDistribution[4]]]

– The defining property for orthogonal matrices:

Inverse[om]===Transpose[om]

– The product of two orthogonal matrices is also orthogonal:

om2=OrthogonalMatrix[RandomVariate[CircularRealMatrixDistribution[4]]]

om.om2

om2.om

– The Kronecker product of two orthogonal matrices is also orthogonal:

KroneckerProduct[om,om2]

KroneckerProduct[om2,om]

Summary

– Structured arrays are an efficient way to represent matrices with inherent structure.

– Use structured arrays instead of lists if the matrices in your application can be represented in terms of these (Vandermonde, Cauchy, Toeplitz, etc.)

THANK YOU!

Questions?

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