Applications of Symbolic Arrays

Itai Seggev, Wolfram Research

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

ABSTRACT

Briefly describe your presentation. By typing your abstract in the box below, you will ensure it is 300 characters or less.

Note: Please Enable Dynamic Content if that button has appeared in the top right of the notebook.

InputField[Dynamic[abstract,
(abstract=If[
StringLength[#]>=300,StringTake[#,300],
#])&],
String,ContinuousAction->True,FieldHint->"Enter your abstract here"]
Text[Dynamic[abstract]]

Plan of attack

Introduce the new functions ArrayExpand and ArraySimplify through survey of matrix and other identities

Show how they can solve problems people report on Stack Exchange

Review symbolic array calculus

Combine the above pieces in 2 detailed examples: least squares and Taylor approximation

Matrix Identities

Many key identities of linear algebra can now be proven and / or applied using ArrayExpand and ArraySimplify (or Simplify/FullSimplify)

The former generally makes things bigger, the latter more compact

Define a couple of n×n matrices:

{a,b}={MatrixSymbol["a",{n,n}],MatrixSymbol["b",{n,n}]}

Basic identities

Definition of inverse:

Simplify[Inverse[a].a==a.Inverse[a]==SymbolicIdentityArray[{n}]]

Expand the inverse of the product:

ArrayExpand[Inverse[a.b]]

Recombine the product:

%//ArraySimplify

Conjugate transpose is involutive:

(a\[ConjugateTranspose])\[ConjugateTranspose]//ArraySimplify

Conjugate transpose expansion:

(a.b)\[Transpose]//ArrayExpand

And back again:

%//ArraySimplify

Inverse and conjugate transpose commute:

ArraySimplify[(Inverse[((ConjugateTranspose[((MatrixSymbol2["a", n, n]))]))])==(ConjugateTranspose[((Inverse[((MatrixSymbol2["a", n, n]))]))])]

Matrices and scalars

Scalar multiplication commutes with trace:

ArrayExpand[Tr[k a],k∈(Complexes)]

Pull a scalar multiple out of a power:

ArrayExpand[MatrixPower[k a,p],k∈(Complexes)]

Determinant of a scalar multiple yields a dimension-dependent multiple:

ArrayExpand[Det[k a],k∈(Complexes)]

Determinant of the inverse is the reciprocal of the determinant:

ArraySimplify[Det[(Inverse[a])]==(1/(Det[a]))]

Determinant of conjugate transpose is the conjugate of the determinant:

really two identities in one!

ArraySimplify[Det[a\[ConjugateTranspose]]==Conjugate[Det[a]]]

Determinant of the matrix exponential is the exponential of the trace:

ArrayExpand[Det[MatrixExp[a]]]

Similarity identities

Similar matrices have the same determinant:

ArrayExpand[Det[b.a.Inverse[b]]]

Do it in the opposite order and with a different command:

ArraySimplify[Det[Inverse[b].a.b]]

Similar matrices have the same trace:

ArraySimplify[Tr[b.a.Inverse[b]]==Tr[(MatrixSymbol2["a", n, n])]]

Matrix exponential generates a one-parameter group:

ArrayExpand[MatrixExp[(s+t)a],(s|t)∈(Reals)]

Vector identities

These functions can deal with vector and general array identities

Will focus on some array identities here; will invoke array identities as needed

Define a couple of vectors:

{v,w}={VectorSymbol["v",{n}],VectorSymbol["w",{n}]}

Some basic identities

Dot products of vectors are commutative:

ArraySimplify[v.w ==w.v]

The definition of transpose (of a matrix / linear operator):

ArraySimplify[v.a.w ==w.a\[Transpose].v]

A vector is its own transpose:

ArraySimplify[v==v\[Transpose]]

A vector is orthogonal to a cross product involving it:

ArrayExpand[v.(v\[Cross]w)]

Note the above made the assumption that n==3; explicit in v×v==(OverVector[0]):

ArraySimplify[v\[Cross]v]

Going a bit fancy

This doesn’t simplify:

ArraySimplify[VectorSymbol["v4",{4}].(VectorSymbol["v4",{4}]\[Cross]VectorSymbol["w4",{4}])]

The correct 4d variant does:

ArraySimplify[(VectorSymbol2["v4", 4]).((VectorSymbol2["v4", 4])\[Cross](VectorSymbol2["w4", 4])\[Cross](VectorSymbol2["v4", 4]))]

Of course, even the term in parentheses is a zero vector:

ArraySimplify[(VectorSymbol2["v4", 4])\[Cross](VectorSymbol2["w4", 4])\[Cross](VectorSymbol2["v4", 4])]

StackExchange

An application is “something that can help someone solve a problem”

I went to the Mathematics Stack Exchange and searched for “matrix algebra” in the linear algebra section

These are the first three problems I found that weren’t either completely abstract or completely concrete

With the new 14.2 features, we can address the problems

Show that 𝐴(((𝐵)^(-1))(𝐴+𝐵)((𝐴)^(-1)))𝐵=𝐴+𝐵

Easy:

ArraySimplify[a.(Inverse[b].(a+b).Inverse[a]).b==a+b]

Show that ((((((Inverse[s]).(SymbolicOnesArray[(ImplicitList[n])]))/((SymbolicOnesArray[(ImplicitList[n])]).(Inverse[s]).(SymbolicOnesArray[(ImplicitList[n])])))))^ᵀ)==(((Inverse[s]).(SymbolicOnesArray[(ImplicitList[n])]))/((SymbolicOnesArray[(ImplicitList[n])]).(Inverse[s]).(SymbolicOnesArray[(ImplicitList[n])]))) for a covariance matrix

The key property of a covariance matrix is that it is symmetric:

s=MatrixSymbol["s",{n,n},Reals,Symmetric[{1,2}]]

ArraySimplify knows what this means:

ArraySimplify[s\[Transpose]]

Then proving the identity is easy:

ArraySimplify[Transpose[(((Inverse[s]).(SymbolicOnesArray[(ImplicitList[n])]))/((SymbolicOnesArray[(ImplicitList[n])]).(Inverse[s]).(SymbolicOnesArray[(ImplicitList[n])])))]==(((Inverse[s]).(SymbolicOnesArray[(ImplicitList[n])]))/((SymbolicOnesArray[(ImplicitList[n])]).(Inverse[s]).(SymbolicOnesArray[(ImplicitList[n])])))]

Show that ∑((((Subscript[v, i])−(OverBar[v])))^2)=vᵀ𝑀v where 𝑀=(Subscript[𝕀, n])-(1/n)(Subscript[1, n])⊗(Subscript[1, n])

Call the re-centered vector z:

z=(v-((v.SymbolicOnesArray[{n}])/n)SymbolicOnesArray[{n}])

Note that (SymbolicOnesArray[(ImplicitList[n])])⊗(SymbolicOnesArray[(ImplicitList[n])]) is simply the matrix of all ones:

FullSimplify[SymbolicOnesArray[{n}]⊗SymbolicOnesArray[{n}]==SymbolicOnesArray[{n,n}]]

ArraySimplify does something, but doesn’t quite get us there:

ArraySimplify[z.z ==v.(SymbolicIdentityArray[{n}]-(1/n)SymbolicOnesArray[{n,n}]).v,n∈(Integers)]

Simplify (and FullSimplify) get us most of the way;

Should be addressed before release

Simplify[z.z ==v.(SymbolicIdentityArray[{n}]-(1/n)SymbolicOnesArray[{n,n}]).v,n∈(Integers)]

We can use ComponentExpand to verify this in any fixed dimension:

Block[{n=10},Expand[ComponentExpand[(((SymbolicOnesArray[(ImplicitList[n,n])]).(VectorSymbol2["v", n]).(VectorSymbol2["v", n])-((Total[(VectorSymbol2["v", n])])^2))/n)]==0]]

For definiteness, ComponentExpand rewrites arrays in terms of components:

Block[{n=2},ComponentExpand[(((SymbolicOnesArray[(ImplicitList[n,n])]).(VectorSymbol2["v", n]).(VectorSymbol2["v", n])-((Total[(VectorSymbol2["v", n])])^2))/n)]]

Calculus review

D can now differentiate with respect to symbolic array variables. For symbolic vectors, this is the familiar gradient, and Grad can be used in the same manner. For higher-rank arrays, this is a generalized gradient. In general, differentiating an array of rank k with respect to an array of rank l will produce an array of rank k+l.

Pre-define a few array variables:

v=VectorSymbol["v",n,Reals];
m=MatrixSymbol["m",{n,n},Reals];
a=ArraySymbol["a",{n,n,n},Reals];

Compute the gradient of a constant, producing a zero vector:

D[5,v]

Do the same computation with Grad:

Grad[5,v]

Compute the derivative of a scalar constant with respect to a matrix:

D[5,m]

The derivative of a constant vector with respect to a matrix; the result has rank three:

D[SymbolicOnesArray[{k}],m]

The derivative of an array variable with respect to itself is always an identity array of appropriate shape:

{D[v,v],D[m,m],D[a,a]}

Differentiate a quadratic form with respect to the vector:

D[v.m.v,v]

Differentiate a vector-valued quadratic form; the result is essentially the same, but expressed differently to align the levels:

D[v.a.v,v]

Differentiate a dot product with respect to the inner array; the result uses a new syntax for Transpose:

D[m.a.v,a]

Confirm that the rank of the result is the sum of the ranks of the two arguments:

TensorRank/@{m.a.v,a,%}

Many array functions have built-in rules:

D[Det[m],m]

And, of course, there is an appropriate chain rule; uses ArrayDot for dotting multiple levels at once:

D[Det[v.a],a]

Higher derivatives are supported:

D[v.v,{v,2}]

As are Div and Laplacian:

Div[v,v]//TensorExpand

Laplacian[v.v,v]//TensorExpand

Least squares

Given a matrix A and vector b, want to find the best solution to Ax==b, i.e., the x that minimizes (Norm[A x-b]).

Define the needed variables:

A=MatrixSymbol["A",{n,m},Reals];
b=VectorSymbol["b",{n},Reals];
x=VectorSymbol["x",{m},Reals];

For a minimum, we need a critical point:

der=D[((Norm[A.x-b])^2),x]

This may not look like the standard expression, but it is:

FullSimplify[der==2A\[Transpose].(A.x-b)]

Solving for x gives x==(Inverse[((Transpose[A]).A)]).A.b==PseudoInverse[A].b assuming (Transpose[A]).A is invertible (for both equal signs)

To guarantee this is a minimum, the second derivative must be positive definite:

D[der,x]

Simplifying, we get the desired invertibility condition:

%//FullSimplify

Taylor Approximation

Use Taylor polynomial approximations of a scalar-valued function of matrices

Expressed in terms of “base point” (MatrixSymbol2["m", n, n]) and “displacement” (MatrixSymbol2["Δm", n, n])

{MatrixSymbol["m",{n,n}],MatrixSymbol["Δm",{n,n}]}

Computing the pieces

Define the function:

Clear[f,m,Δm]
f[m_]:= Det[m]

Compute the first derivative:

D needs to know the shape of the argument, so differentiate with respect to a symbolic matrix and then replace with a variable to define a function

d1[m_]=D[f[(MatrixSymbol2["m", n, n])],(MatrixSymbol2["m", n, n])]/.(MatrixSymbol2["m", n, n])->m

Compute the second derivative:

d2[m_]=D[f[(MatrixSymbol2["m", n, n])],{(MatrixSymbol2["m", n, n]),2}]/. (MatrixSymbol2["m", n, n])->m

This result could be written as follows

tc=(Transpose[(-Det[(MatrixSymbol2["m", n, n])] TensorContract[(Inverse[(MatrixSymbol2["m", n, n])])⊗(Inverse[(MatrixSymbol2["m", n, n])])⊗(SymbolicIdentityArray[(ImplicitList[n,n])]),{{2,5},{3,6}}]+(Inverse[(MatrixSymbol2["m", n, n])])⊗(Transpose[Adjugate[(MatrixSymbol2["m", n, n])]]))])

Verify with ComponentExpand:

Block[{n=4},Expand[ComponentExpand[d2[(MatrixSymbol2["m", n, n])]==tc]]]

Defining the polynomials

The 0th order polynomial is trivial:

t0[m_,Δm_]=f[m]

The first derivative has tensor rank 2; need 2 contractions to get a scalar

In obvious generalization of the Taylor formula from multivariate calculus, use ArrayDot instead of Dot to contract Δm into the derivative

In this particular case, could have written the linear term as TensorContract[d1[m]⊗Δm,{{1,3},{2,4}}]

t1[m_,Δm_]=f[m]+ (ArrayDot[d1[m], Δm, 2])

Similarly, the 2nd derivative has rank 4, and we must dot in two copies of Δm:

Note: even though conceptually this could have been written with TensorContract, it can’t in practice due to the explicitly non-tensorial multiplication in d2. Fixing this friction is on the to-do list.

t2[m_,Δm_]=f[m]+ (ArrayDot[d1[m], Δm, 2]) +(1/2) (ArrayDot[((ArrayDot[d2[m], Δm, 2])), Δm, 2])

Could have been written with the more symmetric form Δm.m.Δm:

Note that equality only holds for square matrices!

t2[(MatrixSymbol2["m", n, n]),(MatrixSymbol2["Δm", n, n])]==f[(MatrixSymbol2["m", n, n])]+ (ArrayDot[d1[(MatrixSymbol2["m", n, n])], (MatrixSymbol2["Δm", n, n]), 2]) +(1/2) (ArrayDot[(MatrixSymbol2["Δm", n, n]), ((ArrayDot[d2[(MatrixSymbol2["m", n, n])], (MatrixSymbol2["Δm", n, n]), 2])), 2])//FullSimplify

Testing the approximation

Generate a matrix and a small perturbation:

a=RandomInteger[{-10^6,10^6},{10,10}]/10^6;
ε=RandomInteger[{-10^6,10^6},{10,10}]/10^8;

This is the numerical value of the exact result

Exact result is a rational with a very large denominator

f[a+ε]//N

This is the 0th order approximation:

Of course identical to evaluating f at a

{t0[a,ε]//N,t0[a,ε]===f[a]}

1st approximation gets a good chunk of the way towards exact result:

t1[a,ε]//N

2nd polynomial gets us 3 or 4 digits

Need to Block n because it literally appears in the result

Normal expands out SymbolicOnesArray and SymbolicZerosArray:

Block[{n=First[TensorDimensions[a]]},Normal[t2[a,ε]]]//N

In conclusion

Version 14.1 gave symbolic arrays (including vectors and matrices) that could be computed with

Version 14.2 allows them to be manipulated symbolically much as version 1 allowed manipulation of elementary algebraic expressions.

A much wider array of problems can now be tackled

Thanks to the rest of the mathematical computation team, especially Adam Strzebonski who did most of the heavy lifting

Questions?