Mike Valenzuela <mickle.mouse at gmail.com> writes:
> [1:multipart/alternative Hide]
> [1/1:text/plain Hide]
> I've been working on a particularly tricky problem and wanted to check my
> math. The problem simply an inverse matrix calculus derivative problem.
> Particularly:
>
> diff( -1/2 *
> transpose(invert(transpose(A)).(z-x)).(invert(transpose(A)).(z-x)), A),
> where A is a matrix ( cholesky decomposition of a covariance matrix), and z
> & x are vectors of the same length.
>
> I ran into several problems:
> 1) Invert only works on matrices if they are given - I need a symbolic
> matrix inversion for any sized matrix. I know differentiation of an
> inverted matrix is possible (see
> orion.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf)
>
> 2) I tried other matrix calculus such as diff (transpose(x) . A . x, x),
> but Maxima returned " 'diff(transpose(x)) . A". The correct answer is " x
> . ( A + transpose(A) ) "
>
> 3) I have tried checking itensor and ctensor packages, but neither of them
> seemed to support a symbolic inverse. I'm not sure that even if it did,
> would I be able to differentiate the resulting expression.
>
> 4) Another thing I looked for, but couldn't find was a symbolic function
> inverse. I know the identity " diff( inverse(f(x)), x) = 1 / ( 'diff(f (
> inverse f( x ) ), x ) " (see
> http://en.wikipedia.org/wiki/Inverse_functions_and_differentiation)
>
>
> So I guess this is a feature request unless someone can tell me how to use
> maxima to differentiate an expression (involving matrix inverses) with
> respect to the matrix. Features specifically useful would be:
> a) "Full" symbolic matrix inversion, which operates on any sized matrix and
> does not require specifying its components.
>
> b) Matrix calculus
>
> c) Symbolic function inverses (along with the associated calculus rules)
>
> or
>
> d) Better itensor/centsor documentation demonstrating how to handle a
> problem similar to mine.
Maxima has the capability to have simplification rules added to its
built-in simplifier. Here are some:
(%i1) kill(all)$
(%i1) display2d:false $
(%i2) matchdeclare(X,lambda([u],get(u,matrixp))) $
(%i3) matchdeclare(t,atom) $
(%i4) put(B,true,matrixp) $
(%i5) depends([A,B],x) $
(%i6) tellsimp('diff(X^^(-1),t),-X^^(-1) . diff(X,t) . X^^(-1));
(%i7) tellsimp('diff('transpose(X),t),transpose(diff(X,t)));
tellsimp is a command that adds simplification rules; matchdeclare lets
us apply these rules selectively.
Now, we can try things out:
(%i8) diff(B,x);
(%o8) 'diff(B,x,1)
(%i9) diff(B^^(-1),x);
(%o9) -B^^(-1) . ('diff(B,x,1)) . B^^(-1)
The correct answer.
(%i10) diff(A,x);
(%o10) 'diff(A,x,1)
(%i11) diff(A^^(-1),x);
(%o11) 'diff(A^^(-1),x,1)
Because we have not told Maxima that A is a matrix, the simplifier
leaves the expression unevaluated.
(%i12) diff(transpose(B),x);
(%o12) 'transpose('diff(B,x,1))
Again, the correct answer.
(%i13) diff(transpose(B). A . B, x);
(%o13) 'transpose('diff(B,x,1)) . A . B+'transpose(B) . ('diff(A,x,1))
. B
+'transpose(B) . A . 'diff(B,x,1)
Looks good, to..
(%i14) diff(transpose(B). C . B, x);
(%o14) 'transpose('diff(B,x,1)) . C . B+'transpose(B) . C . 'diff(B,x,1)
Looks good (Maxima treats C as a constant, unlike in %o13).
--
Leo Butler <l_butler at users.sourceforge.net>
SDF Public Access UNIX System - http://sdf.lonestar.org