Oh and also a thank you to Robert (sorry, I missed your name in my original
thank you list, it was an accident)
On Thu, Feb 23, 2012 at 1:14 PM, Mike Valenzuela <mickle.mouse at gmail.com>wrote:
> 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
>>>
>>
>>
>