Menu
▾
▴
Re: [Maxima-discuss] How to get a constant factor 'out' of a matrix?
|
From: Stavros M. <mac...@gm...> - 2023-07-10 22:57:29
|
How about this? The test function *testfmd *should return all zeros if the
decomposition is correct. The decomposition is (unsurprisingly) not unique.
In some cases, my code gives different results from your code. In
particular, yours can give fractions in the middle matrix, e.g.,
for factormatrixgcds(matrix([2,2],[2,1])).
-s
vecgcd(l):=block([g:lreduce('gcd,l),listarith:true],
if g=0 then [0,0*l] else [g,l/g])$
factormatrixdiags(m):=
block( [rowgcd,rows,colgcd,cols,rowdiag,coldiag,rowdiaggcd,coldiaggcd],
if not matrixp(m) then error("factormatrixdiags requires a matrix
argument:",m),
rowgcd: map('vecgcd, args(m)),
rows: funmake('matrix,map('second,rowgcd)),
colgcd: map('vecgcd,args(transpose(rows))),
cols: funmake('matrix,map('second,colgcd)),
rowdiag: map('first,rowgcd),
coldiag: map('first,colgcd),
rowdiaggcd: vecgcd(rowdiag),
coldiaggcd: vecgcd(coldiag),
[first(rowdiaggcd)*first(coldiaggcd)/max(1,gcd(first(rowdiaggcd),first(coldiaggcd))),
apply('diag_matrix,second(rowdiaggcd)),
transpose(cols),
apply('diag_matrix,second(coldiaggcd))])$
(%i1) load("~/maxima/mgcdx.mac");
(%o1) ~/maxima/mgcdx.mac
(%i2)
testfmd(m):=(print(""),disp(m),disp(gg:factormatrixdiags(m)),apply(".",gg)-m)$
(%i3) testfmd(matrix([0,0],[0,0]));
[ 0 0 ]
[ ]
[ 0 0 ]
[ 0 0 ] [ 0 0 ] [ 0 0 ]
[0, [ ], [ ], [ ]]
[ 0 0 ] [ 0 0 ] [ 0 0 ]
[ 0 0 ]
(%o3) [ ]
[ 0 0 ]
(%i4) testfmd(matrix([0,0],[1,0]));
[ 0 0 ]
[ ]
[ 1 0 ]
[ 0 0 ] [ 0 0 ] [ 1 0 ]
[1, [ ], [ ], [ ]]
[ 0 1 ] [ 1 0 ] [ 0 0 ]
[ 0 0 ]
(%o4) [ ]
[ 0 0 ]
(%i5) testfmd(matrix([5,0],[0,5]));
[ 5 0 ]
[ ]
[ 0 5 ]
[ 1 0 ] [ 1 0 ] [ 1 0 ]
[5, [ ], [ ], [ ]]
[ 0 1 ] [ 0 1 ] [ 0 1 ]
[ 0 0 ]
(%o5) [ ]
[ 0 0 ]
(%i6) testfmd(matrix([72,60],[630,336]));
[ 72 60 ]
[ ]
[ 630 336 ]
[ 2 0 ] [ 2 5 ] [ 3 0 ]
[6, [ ], [ ], [ ]]
[ 0 7 ] [ 5 8 ] [ 0 1 ]
[ 0 0 ]
(%o6) [ ]
[ 0 0 ]
(%i7) testfmd(matrix([72,60],[0,336]));
[ 72 60 ]
[ ]
[ 0 336 ]
[ 1 0 ] [ 1 5 ] [ 6 0 ]
[12, [ ], [ ], [ ]]
[ 0 28 ] [ 0 1 ] [ 0 1 ]
[ 0 0 ]
(%o7) [ ]
[ 0 0 ]
(%i8) testfmd(matrix([72,60],[630,336])^^3);
[ 18517248 10804320 ]
[ ]
[ 113445360 66056256 ]
[ 2 0 ] [ 7144 12505 ] [ 3 0 ]
[432, [ ], [ ], [ ]]
[ 0 7 ] [ 12505 21844 ] [ 0 1 ]
[ 0 0 ]
(%o8) [ ]
[ 0 0 ]
(%i15) testfmd(matrix([6,120,9],[8,16,4]));
[ 6 120 9 ]
[ ]
[ 8 16 4 ]
[ 2 0 0 ]
[ 3 0 ] [ 1 10 3 ] [ ]
[1, [ ], [ ], [ 0 4 0 ]]
[ 0 4 ] [ 1 1 1 ] [ ]
[ 0 0 1 ]
[ 0 0 0 ]
(%o15) [ ]
[ 0 0 0 ]
(%i16) testfmd(transpose(matrix([6,120,9],[8,16,4])));
[ 6 8 ]
[ ]
[ 120 16 ]
[ ]
[ 9 4 ]
[ 2 0 0 ] [ 1 2 ]
[ ] [ ] [ 3 0 ]
[1, [ 0 8 0 ], [ 5 1 ], [ ]]
[ ] [ ] [ 0 2 ]
[ 0 0 1 ] [ 3 2 ]
[ 0 0 ]
[ ]
(%o16) [ 0 0 ]
[ ]
[ 0 0 ]
(%i17) testfmd(matrix([72,60],[0,336])^^-1);
[ 1 5 ]
[ -- - ---- ]
[ 72 2016 ]
[ ]
[ 1 ]
[ 0 --- ]
[ 336 ]
1 [ 1 0 ] [ 1 - 5 ] [ 28 0 ]
[----, [ ], [ ], [ ]]
2016 [ 0 6 ] [ 0 1 ] [ 0 1 ]
[ 0 0 ]
(%o17) [ ]
[ 0 0 ]
On Sun, Jul 9, 2023 at 1:31 PM Henry Baker <hb...@pi...> wrote:
> Oops! Still issues with division by zero.
>
>
>
> Hopefully attached is the last&best version.
>
>
>
> -----Original Message-----
> From: Henry Baker <hb...@pi...>
> Sent: Jul 9, 2023 8:49 AM
> To: Stavros Macrakis <mac...@gm...>
> Cc: <max...@li...>
> Subject: Re: [Maxima-discuss] How to get a constant factor 'out' of a
> matrix?
>
>
>
> Whoops!
>
>
>
> That code didn't work properly in the case of an overall factor in the
> matrix M.
>
>
>
> Attached is the fixed code for
>
>
>
> /* factormatrixgcds(M) returns a 4-element list: g, D, N, E */
>
>
> /* g is overall gcd of all elements of matrix M,
> D is diagonal matrix of row gcds of M,
> E is diagonal matrix of col gcds of M,
> N is the result after removing g, row gcds, and col gcds from M.
>
> M = g*(D.N.E)
>
> */
>
>
>
> -----Original Message-----
> From: Henry Baker <hb...@pi...>
> Sent: Jul 9, 2023 7:16 AM
> To: Stavros Macrakis <mac...@gm...>
> Cc: <max...@li...>
> Subject: Re: [Maxima-discuss] How to get a constant factor 'out' of a
> matrix?
>
>
>
> Re: why don't you submit your code?
>
>
>
> Here it is (I apologize for the bad formatting); do with it whatever you
> like:
>
>
>
> mtest:transpose([a,b,c,d]).[e,f,g];
>
>
>
>
> factormatrixgcds(M):=block([nrows,ncols,i,j,rowgcds,rowgcd,colgcds,colgcd,N],
> [nrows,ncols]:matrix_size(M),
> rowgcds:[],
> for i:1 thru nrows do
> (rowgcd:0,
> for j:1 thru ncols do
> rowgcd:gcd(rowgcd,M[i,j]),
> push(rowgcd,rowgcds)),
> rowgcds:reverse(rowgcds),
> colgcds:[],
> for j:1 thru ncols do
> (colgcd:0,
> for i:1 thru nrows do
> colgcd:gcd(colgcd,M[i,j]),
> push(colgcd,colgcds)),
> colgcds:reverse(colgcds),
> N:copymatrix(M),
> for i:1 thru nrows do
> for j:1 thru ncols do
> N[i,j]:N[i,j]/(rowgcds[i]*colgcds[j]),
>
> [apply('diag_matrix,rowgcds),N,apply('diag_matrix,colgcds)]);
>
>
>
> (%i1) mtest;
> [ a e a f a g ]
> [ ]
> [ b e b f b g ]
> (%o1) [ ]
> [ c e c f c g ]
> [ ]
> [ d e d f d g ]
> (%i2) factormatrixgcds(mtest);
> [ a 0 0 0 ] [ 1 1 1 ]
> [ ] [ ] [ e 0 0 ]
> [ 0 b 0 0 ] [ 1 1 1 ] [ ]
> (%o2) [[ ], [ ], [ 0 f 0 ]]
> [ 0 0 c 0 ] [ 1 1 1 ] [ ]
> [ ] [ ] [ 0 0 g ]
> [ 0 0 0 d ] [ 1 1 1 ]
> (%i3) lreduce(".",%);
> [ a e a f a g ]
> [ ]
> [ b e b f b g ]
> (%o3) [ ]
> [ c e c f c g ]
> [ ]
> [ d e d f d g ]
>
>
>
> -----Original Message-----
> From: Stavros Macrakis <mac...@gm...>
> Sent: Jul 7, 2023 9:20 AM
> To: <hb...@pi...>
> Cc: Barton Willis <wi...@un...>, <max...@li...
> >
> Subject: Re: [Maxima-discuss] How to get a constant factor 'out' of a
> matrix?
>
>
> Why don't you submit your code to *contrib*?
>
> On Fri, Jul 7, 2023 at 11:36 AM Henry Baker <hb...@pi...> wrote:
>
>> In the past, I have written a bit of code that *factors* any matrix M
>>
>> into a product M = D.N.E, where D,E are *diagonal* matrices whose
>>
>> elements are the corresponding gcd's of the rows/columns of M.
>>
>>
>>
>> Note that this factorization a) is totally rational; and b) not very
>> expensive
>>
>> to calculate.
>>
>>
>>
>> In some of my calculations, the D & E matrices have a reasonable amount
>>
>> of 'stuff' that isn't very important to understanding what's really going
>> on,
>>
>> so I can focus all my attention on the 'N' matrix that's left after these
>>
>> irrelevant factors are removed.
>>
>>
>>
>> I'd love to see such a matrix factorization as a standard function in
>> Maxima.
>>
>>
>>
>> -----Original Message-----
>> From: Stavros Macrakis <mac...@gm...>
>> Sent: Jul 7, 2023 6:26 AM
>> To: Barton Willis <wi...@un...>
>> Cc: max...@li... <
>> max...@li...>
>> Subject: Re: [Maxima-discuss] How to get a constant factor 'out' of a
>> matrix?
>>
>>
>> In your case, as Barton pointed out, there are options to keep the
>> determinant as a multiplier of the matrix.
>>
>> More generally, if you have an arbitrary matrix with a constant factor
>> that you want to pull out, you can do something like this:
>>
>> (%i6) mm:matrix([2*x^2*y, 4*x], [8*x, 16*x*y]);
>>
>> [ 2 ]
>> (%o6) [ 2 x y 4 x ]
>> [ ]
>> [ 8 x 16 x y ]
>>
>> (%i7) apply('append,args(mm)) /* append all the rows into one list */ ;
>>
>> 2
>> (%o7) [2 x y, 4 x, 8 x, 16 x y]
>>
>> (%i8) xreduce('gcd,%) /* IMPORTANT to quote 'gcd */ ;
>>
>> (%o8) 2 x
>>
>> (%i9) mm/%;
>>
>> [ x y 2 ]
>> (%o9) [ ]
>> [ 4 8 y ]
>>
>> On Fri, Jul 7, 2023 at 6:48 AM Barton Willis via Maxima-discuss <
>> max...@li...> wrote:
>>
>>> Hi Wolfgang,
>>>
>>> Try setting the option variables * detout* to true*, doallmxops* to
>>> false, and *doscmxop*s to false. The user documentation for invert
>>> (input: ? invert) mentions this
>>>
>>> * When 'detout' is 'true', the determinant is factored out of the*
>>> * inverse. The global flags 'doallmxops' and 'doscmxops' must be*
>>> * 'false' to prevent the determinant from being absorbed into the*
>>>
>>>
>>> Example
>>>
>>> (%i19) (detout : false, doallmxops : false, doscmxops : false);
>>> (%o19) false
>>>
>>> (%i20) invert(matrix([a_11,-a_12], [-a_12,a_22]));
>>> (%o20) matrix([a_22, a_12],[a_12, a_11])/(a_11*a_22-a_12^2)
>>>
>>>
>>> If you do this frequently, you might like to define your own
>>> function—something like
>>>
>>> (%i1) my_invert(M) := block([detout : true,doallmxops : false, doscmxops
>>> : false], invert(M))$
>>>
>>> (%i2) my_invert(matrix([42,4],[8,107]));
>>> (%o2) matrix([107, -4], [-8, 42])/4462
>>>
>>> Let us know if this doesn't resolve your question or if you have other
>>> questions.
>>>
>>> --Barton
>>> ------------------------------
>>> *From:* Wolfgang Hugemann <Au...@Hu...>
>>> *Sent:* Friday, July 7, 2023 03:26
>>> *To:* max...@li... <
>>> max...@li...>
>>> *Subject:* [Maxima-discuss] How to get a constant factor 'out' of a
>>> matrix?
>>>
>>> Non-NU Email
>>>
>>> Dear list members,
>>>
>>> after years of struggling with Maxima on my own, I will give this
>>> mailing list a try. In general, Maxima does an amazing job, but I often
>>> struggle with operations that seem to be most simple.
>>>
>>> These two lines
>>>
>>> matrix([a_11,-a_12], [-a_12,a_22]);
>>> M:invert(%);
>>>
>>> will yield
>>>
>>> (%o2)
>>>
>>> matrix([a_22/(a_11*a_22-a_12^2),a_12/(a_11*a_22-a_12^2)],[a_12/(a_11*a_22-a_12^2),a_11/(a_11*a_22-a_12^2)])
>>>
>>> How do I get the common constant factor (a_11*a_22-a_12^2) 'out' of the
>>> matrix?
>>>
>>> Wolfgang Hugemann
>>>
>>>
>>>
>>> _______________________________________________
>>> Maxima-discuss mailing list
>>> Max...@li...
>>>
>>> https://urldefense.com/v3/__https://lists.sourceforge.net/lists/listinfo/maxima-discuss__;!!PvXuogZ4sRB2p-tU!DRR2pH5_PAf9nfmgsPbH0p567ITTFL4MUdQYN1VOoonNiMz0BktY5FFh5JZWUMkhZ9qwZ8AOTY6kY8E$
>>> _______________________________________________
>>> Maxima-discuss mailing list
>>> Max...@li...
>>> https://lists.sourceforge.net/lists/listinfo/maxima-discuss
>>>
>>
>>
>
>
>
>
>
>
|
