So I decided to start setting up a drop. For now this box is tied to my
email address (but I suppose we can change at any point in time, I really
don't care and don't use drop box for anything else yet). I can send an
invite to share the project folder with anyone that is interested.
On Thu, Feb 23, 2012 at 1:20 PM, Mike Valenzuela <mickle.mouse at gmail.com>wrote:
> 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
>>>>
>>>
>>>
>>
>