Help building a matrix-calculus package



On 2/22/12, Mike Valenzuela <mickle.mouse at gmail.com> wrote:

> (A_1 . A_2 ... A_n)^^-1 = A_n^^-1 ... A_2^^-1 . A_1^^-1

matchdeclare (aa, nctimesp);
nctimesp (e) := not atom(e) and op(e) = ".";
matchdeclare (xx, all);
tellsimp (aa^^xx, map (lambda ([e], e^^xx), aa));

(a . b . c)^^-1;
 => a^^(-1) . b^^(-1) . c^^(-1)

> Also, there are some recursive rules (a recursive rule for just 2
> elements could work for the above problem), which may work well using
> tellsimp, but I don't know.

At present, Maxima's pattern matcher can't deduce rules for
n-ary operators from a pattern with just two operands of the same type.
It would be useful and interesting to do that or other kinds of
deductions based on properties of the operator (commutative, etc).
It would be some work (I can do it) but not really a big job, so if you're
interested, let's talk it over.

> partial( g(U) ) / partial(X_ij) = Trace( transpose(partial g(U) /
> partial U) . partial U / partial X_ij),

Dunno, didn't try this one yet.

> partial ( transpose(X) . B . X ) / partial (X_ij) = transpose(X) . B .
> single_entry_matrix(i,j) + single_entry_matrix(j,i) . B . X,
> where the single_entry_matrix size is determined by the expression (in
> this case, it is the same size as transpose(X) . B . X)

Let E = single_entry_matrix. (I didn't try to write this function.)

matchdeclare ([xx, bb], nonscalarp);
matchdeclare ([ii, jj], all);
/* let built-in simplifier act on diff(transpose(xx) ...),
 * then construct new rule for simplified expression
 */
'diff (transpose (xx) . bb . xx, xx [ii, jj]);
tellsimpafter (''%, transpose(xx) . bb . E(ii, jj) + E(jj, ii) . bb . xx);

declare ([X, B], nonscalar);
diff (transpose (X) . B . X, X[i, j]);
 => 0
'diff (transpose (X) . B . X, X[i, j]);
 =>  transpose(X) . B . E(i, j) + E(j, i) . B . X

There are various confusing details here (sorry).
Would it work as well without the nonscalarp test? Maybe.
Why is diff(...) different from 'diff(...) ? The rule was written for the
noun expression 'diff(...), it doesn't match diff(...); maybe it should.

> Is there some way to define a variable as a row or column vector? It
> would be useful for expressions like:
> partial transpose( b ) . transpose(X) . X . c / partial X = X( b .
> transpose(c) + c . transpose(b) ),
> specifically I need to know that b and c are column vectors for the
> above identity to work.

Well, a conventional way would be to represent them as matrices
with 1 row or 1 column. That's not entirely satisfactory from an
algebraic point of view ... I don't like conflating different types of
objects for mere convenience.

A more principled approach is to define a vector type and make it
possible to have symbolic vectors (e.g. declare(x, vector)) in addition
to literal vectors (e.g. x : vector(a, b, c)). I don't know if someone has
already done so; several matrix/vector packages have been written
over the years, maybe one of them has that stuff.

> How can I tell Maxima to "Rotate" matrices inside a trace?
> Tr( A . B . C . single_entry_matrix(i, j) . D . E) = Tr( D . E . A . B
> . C . single_entry_matrix(i, j) ) = the ith, jth index of ( D . E . A
> . B . C )

I don't see an obvious way to do it. Certainly it can be done in a
messy way, but I don't see a neat way to do it.

> Thanks in advance for any advice. I know this is a very ambitious
> undertaking since the first two chapters contain about 100 identities.
> I suspect at least some can be re-derived using the chain rule and
> some of the basic rules. Again, I am not sure how best to go about
> building a package to really help me out.

Well, I think the way to go about this is to write out the list of
identities you want, then implement as many of them as you can,
maybe reusing bits from existing share packages and maybe
punting to built-in rules sometimes. (Note that I'm specifically
recommending that you make your list first and review existing
packages afterwards.)

I would be happy to help write rules and I'm willing to consider
extensions to the pattern matching code.

all the best

Robert Dodier