Thank you Stavros, Richard, and Barton,
I think I may create either a Google doc or a drop-box project enumerating
the rules I wish to implement, a link to the publicly available book
"Matrix Cookbook (2008 ed.)", and a code file (dropbox has a version
control system). I will then make it freely available to everyone who would
want to help contribute. I typically have more time available for this at
night (but I might be able to convince my academic adviser that this
project would aid my research), so this may not be up for another day or so.
Also, since I am a graduate student, I believe I might as well poll my
university's math department to see if anyone there knows about Maxima and
would be interested in helping (they may be able to publish a paper
discussing the package). I believe it would be invaluable to have one or
more experts generate test-cases / test-problems.
One side note - I believe it ought to be possible to enter the matix
calculus chain rule, even if it is not entered recursively, because Maxima
appears to have a chain rule for traditional calculus. I don't know how it
is implemented, but the matrix calculus chain rule is very similar and
should be able to reuse much of the same logic.
Once again, thank you all.
On Thu, Feb 23, 2012 at 10:56 AM, Stavros Macrakis <macrakis at alum.mit.edu>wrote:
> Mike,
>
> Robert's general approach works, but the rule itself has a bug.
>
> I would suggest something like this:
>
> matchdeclare (ncprod, nctimesp);
> nctimesp (e) := not atom(e) and op(e) = ".";
> matchdeclare(negint,lambda([x],featurep(x,integer) and x<0));
> tellsimp (ncprod^^negint, (map (lambda ([e], e^^-1), reverse(aa) ) ^^
> -negint));
>
> The next issue you will encounter is simplifying expressions like
>
> (a.b.b)^^2 . (b^^-1 . a^^-1)^^5 . a.b.a
>
> Here's a simple way to do that:
>
> ncexpandsimp(ex):=
> block([dotexptsimp:true],
> expand(block([dotexptsimp:false],expand(ex))))$
>
> ncexpandsimp( (a.b.b)^^2 . (b^^-1 . a^^-1)^^5 . a.b.a ) =>
> a . b^^2 . a . b . (a^^(-1) . b^^(-1))^^3
>
> Of course, there is more than one way to write some expressions, so for
> example
>
> ncexpandsimp((a.b.a)^^2 . b) => a.b.a.(a.b)^^2
>
> which may not be what you want.
>
> If you want to handle symbolic exponents, e.g. (a.b.a)^^k . a^^-1 =>
> (a.b.a)^^(k-1).a.b, you'll need to write additional simplification rules.
>
> -s
>
>
> On Thu, Feb 23, 2012 at 10:57, Robert Dodier <robert.dodier at gmail.com>wrote:
>
>> 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
>> _______________________________________________
>> Maxima mailing list
>> Maxima at math.utexas.edu
>> http://www.math.utexas.edu/mailman/listinfo/maxima
>>
>
>