Calculus and Algebra with Symbolic Arrays

Adam Strzeboński, Wolfram Research

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

ABSTRACT

Mathematica 14.1 introduces new functionality for computation with symbolic arrays. It includes symbolic representation of arrays, computation of derivatives of array functions with respect to symbolic array variables, as well as support for array variables in several symbolic computation functions. The upcoming release of Mathematica will add functionality for transforming and simplifying symbolic array expressions and extend support for array variables to more functions. In my talk, I will discuss the new and upcoming symbolic array functionality and show some examples.

Symbolic Arrays

NonThreadable attribute

New attribute that prevents threading of listable functions. It is given to symbols that may be used to construct expressions that represent arrays.

First/@Select[ToExpression[#,InputForm,Hold]&/@Names["System`*"],MemberQ[Attributes@@#,NonThreadable]&]

Symbolic array variables

• VectorSymbol

• MatrixSymbol

• ArraySymbol

Special arrays

• SymbolicZerosArray

• SymbolicOnesArray

• SymbolicIdentityArray

• SymbolicDeltaProductArray

Symbolic array calculus

D works with symbolic array functions and variables.

Symbolic array transformations

• ArrayExpand

• ArraySimplify

• ComponentExpand

• TensorExpand

• TensorReduce

Symbolic Array Variables and Threading

A long-standing difficulty with using symbolic vector, matrix, and array variables was that evaluation of arithmetic operations treated symbolic expressions as scalars and threaded them into list arguments.

x+{1,2}

For instance, we would like to set up a convex optimization constraint a.x+b\[VectorGreaterEqual]0. Because of threading we could not do it directly.

a={{1,2},{1,0}};
b={-3,1};
a.x+b\[VectorGreaterEqual]0

As suggested in the documentation of ConvexOptimization, we needed to use Inactive[Plus].

Inactive[Plus][a.x,b]\[VectorGreaterEqual]0

ConvexOptimization[{1,1}.x,%,x]

NonThreadable attribute

Mathematica 14.1 introduced a new NonThreadable attribute which specifies that a symbol, and expressions that have this symbol as the head, should not be threaded into lists.

SetAttributes[sym,NonThreadable]

sym+{1,2}

sym[1,2,3]+{1,2,3}

Most standard built-in functions that produce array results now have the NonThreadable attribute.

Inverse[mat]+{{1,2},{3,4}}

Symbolic Array Variables

VectorSymbol, MatrixSymbol, and ArraySymbol are new objects that can be used for symbolic representation of vectors, matrices, and arrays.

v=VectorSymbol["v",2]

The new objects have the NonThreadable attribute.

v+{1,2}

Now we can set up the convex optimization constraint directly.

a.v+b\[VectorGreaterEqual]0

ConvexOptimization[{1,1}.v,%,v]

Array variables can be used in solvers.

a=MatrixSymbol["a",{2,2}]

Solve[a.{1,2}=={3,4}&&Det[a]==5&&Tr[a]==6,a]

a=ArraySymbol["a",{2,2,2}]

Minimize[{Tr[a],a^2<=2},a]

Symbolic Array Expressions

Represent a matrix.

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

Construct a symbolic array expression.

TensorProduct[a.Transpose[a],ConjugateTranspose[a].a]

TensorDimensions[%]

Represent an array of depth 3.

a=ArraySymbol["a",{p,q,r}]

a-a gives a SymbolicZerosArray.

a-a

Adding one results in SymbolicOnesArray.

%+1

With numeric dimensions Normal gives an explicit array.

Normal[%/.{p->2,q->3,r->4}]

Array Derivatives

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

Derivative of a vector valued function wrt. a scalar variable.

D[a[x].v[x],x]

Derivative of a vector valued function wrt. a vector variable.

D[a.v,v]

Derivative of a vector valued function wrt. a matrix variable.

D[a.v,a]

TensorDimensions[%]

SymbolicIdentityArray[dims] is the result of differentiating an array variable with dimension dims wrt. itself.

D[a,a]

TensorDimensions[%]

Derivative of a matrix valued function wrt. a matrix variable.

D[Inverse[a],a]

Specify the domain of a vector variable.

u=VectorSymbol["u",n,Reals];

Variance is differentiable only for real vectors.

D[Variance[u],u]

Chain Rule and ArrayDot

f=MatrixSymbol["f",{m,n}];
g=ArraySymbol["g",{p,q,r}];
v=VectorSymbol["v",k];

Derivative of a composition of symbolic array-valued functions wrt. an array variable.

d=D[f[g[v]],v]

{Head[d],List@@d}

ArrayDot[a,b,k] computes the product of arrays a and b obtained by summing up products of terms over the last k dimensions of a and the first k dimensions of b.

TensorDimensions[d[[1]]]

TensorDimensions[d[[2]]]

TensorDimensions[d]

Example: Taylor Approximations

Find the first three terms of the Taylor series of Inverse[a] at c=(Subscript[I, n]).

a=MatrixSymbol["a",{n,n}];
f=Inverse[a];
c=SymbolicIdentityArray[{n}];

The first term is f(c).

f/.a->c

Simplify the result with ArraySimplify (coming in the next version).

t0=ArraySimplify[%]

The second term is f'(c)(Underscript[·, 2])(a-c).

D[f,a]/.a->c

t1=ArrayDot[ArraySimplify[%],a-c,2]

The third term is (1/2)f''(c)(Underscript[·, 2])(a-c)(Underscript[·, 2])(a-c).

D[f,{a,2}]/.a->c

t2=1/2ArrayDot[ArrayDot[ArraySimplify[%],a-c,2],a-c,2]

Compare order 0, 1, and 2 Taylor approximations for n=7 and a with random entries within (10^(-10)) of the entries of the identity matrix.

aval=IdentityMatrix[7]+RandomReal[{-10^-10,10^-10},{7,7},WorkingPrecision->50];

Norm[Normal[f-t0/.{a->aval,n->7}]]

Norm[Normal[f-t0-t1/.{a->aval,n->7}]]

Norm[Normal[f-t0-t1-t2/.{a->aval,n->7}]]

ArrayExpand

The upcoming version of Mathematica will contain a new function ArrayExpand which expands out symbolic array operations. The function makes use of multi-linearity of array operations, as well as of numerous array, matrix, and vector operation identities.

Clear[{a,b,c,d,e,f,u,v}]

ArrayExpand[(a+b).(c+d)]

ArrayExpand[Tr[s a+t b],Element[s|t,Complexes]]

ArrayExpand[(Inverse[((MatrixSymbol2["a", n, n]).(MatrixSymbol2["b", n, n]))])]

ArrayExpand[Transpose[a.b.c],Element[a|b|c,Matrices[{n,n}]]]

ArrayExpand[Det[3(MatrixSymbol2["a", n, n])]]

ArrayExpand[Cross[a,b,c].Cross[d,e,f]]

Example: Formula Derivation

Derive a formula for the gradient of the least-squares cost function (example from Mathematica StackExchange post).

X=MatrixSymbol["X",{m,n}];
y=VectorSymbol["y",m];
θ=VectorSymbol["θ",n];

Define the cost function.

J=1/2(X.θ-y).(X.θ-y)

Compute the gradient.

g=Grad[J,θ]

Simplify the gradient.

ArrayExpand[g]

ComponentExpand

ComponentExpand expands out symbolic array expressions in terms of indexed components.

ComponentExpand[VectorSymbol["v",3]]

ComponentExpand[Det[a.b],Element[a|b,Matrices[{2,2}]]]

ComponentExpand[VectorLess[{u,v}],Element[u|v,Vectors[3]]]

ComponentExpand[SymbolicOnesArray[{2,2}].MatrixSymbol["a",{2,2}]]

Thank you for listening!