Harry Garst

279 Reputation

6 Badges

19 years, 331 days

MaplePrimes Activity


These are replies submitted by Harry Garst

@dharr
It seems to work. You’ve done some great work!
But even when using random matrices, it is very slow. I am afraid this method creates a lot of redundancies in the calculation. The selection matrices involved only affect the inner matrix. Maybe it is similar to what I did before for a single matrix.

Inverse_single_matrix.mw

@dharr 

That is very interesting. Tomorrow I will check whether it is generalizable.
Thanks a lot for your help.

Harry

@dharr 

Indeed, I also found that generalizing the three matrices to other dimensions is much more difficult. Unfortunately, I do not understand why it only works in special cases. However, I did find a generalization for the case where A is 4×5, B is 5×5, and C is 4×5. The formula is an extension of the previous one, but not a complete one. The solution uses only the first 100 of the 1000 cases in the sequence; the last 1, 2, and 3 remain unchanged.

My interest in the inverse of a product of a rectangular matrix, an invertible square matrix, and another rectangular matrix comes from statistics: for calculating standard errors, the asymptotic covariance matrix has to be inverted. In this setting, the middle matrix may be very large, but its inverse is very easy to obtain. However, the presence of the two rectangular matrices makes the overall inversion very difficult symbolically (numerically it is, of course, not a problem).
Now that I am retired I can play with these problems without worrying about its scientific value.

A4_B5_C4.mw

@dharr 

I think this is the Cauchy–Binet approach. It also breaks the problem into small pieces, but it looks quite different at first sight.
Cauchy_Binet.mw

@dharr 

restart; with(LinearAlgebra)

dim := 6; bim := 36

A := RandomMatrix(bim, bim, shape = symmetric)

B := RandomMatrix(dim, bim)

X := B.A.B^%T

slim := bim-dim

C := `<,>`(`<,>`(B), `<|>`(ZeroMatrix(slim, dim), IdentityMatrix(slim)))

Y := C.A.C^%T

simplify(SubMatrix(Y, 1 .. dim, 1 .. dim)-X)

Matrix(%id = 36893489527472642156)

(1)

simplify(SubMatrix(eval(1/Y), 1 .. dim, 1 .. dim)-SubMatrix(eval(1/Y), 1 .. dim, dim+1 .. bim).(1/SubMatrix(eval(1/Y), dim+1 .. bim, dim+1 .. bim)).SubMatrix(eval(1/Y), dim+1 .. bim, 1 .. dim)-1/X)

Matrix(%id = 36893489527538739188)

(2)


Download Augment_rectangular_matrix_to_square_matrix.mw

@dharr 

Thank a lot! Using your formula it is easy to extend the model to 6 and 8 dimensions of the inner square matrix.
Adjoint383_3Dharr.mw

Another way of inverting a matrix product including rectangular matrices is to augment these matrices into invertable matrices and then use the Block Inversion formula in reverse.
See next post

Thanks for all your answers! Now, I can finish this project.
Harry Garst

initializing the array as follows:

visited := Array(1 .. n, fill = false)

and now it works!!

I think I like Gemini!!! (and it is very polite)

Gemini:

You are absolutely correct! Using fill = false is the more robust and recommended way to initialize a Maple Array with a specific value. While Array(1..n, false) sometimes works due to Maple's implicit type coercion, fill = false is the explicitly correct and consistent method.

Here's the corrected and improved PermutationParity procedure:

@dharr 

Thanks for your suggestion. However, if problems are too complex, I also fall back on numerical methods.

You wrote: "The main issue here is that the inverse of a matrix with symbolic entries is too complicated to be useful, except for very small matrices."

In my experience, it is possible to invert and interpret rather large and complicated matrices. All depends on the extent the matrices are structured or patterned, which they often are in statistics. And of course it depends on the power of the computer you are working on (my PC has 64GB of internal memory and a rather fast CPU). I can give you some spectacular results of Maple's capabilities in inverting and factoring large complicated matrices (of course with highly structured). But if the previous Kronecker example needs too much time, it will not work on your computer.
The main topic of my question has not been answered. Can someone show me how to include basis math results into the symbolic calculation in Maple. After all, no computation was needed if basic math knowledge was utilized.
Maybe as assumptions. Or using AI?

@sursumCorda 

Thanks a lot! You are the winner! And I learned something new: $ command!!

But for a matrix of size 8 the following is even faster than the build in command for determinant:
test_Determinants_symbolic_Modified_8_dimensions_Modified.mw

It always happens to me if I have too much new output in my worksheet. Deleting large parts of output solves the problem. 

@Carl Love 

Thanks a lot. Without the help of the experts of Mapleprimes it would be difficult to figure this out. You may say RTFM, but the Maple Help System is not always very helpful. Sometimes I suspect that a paragraph in the Maple Help System is only comprehensible for someone who is already an expert on the subject.
Maybe after clarifying a user question on Mapleprimes, a small addendum where the user could find all the information him- or herself, would be convenient. :-). If not present, update the help file.

@acer 

Thanks for your explanation. I still do not understand the following result:

``

restart; with(LinearAlgebra)

komb := proc (n::integer, m::integer)::indexed; local a, L, j; L := [seq(a[j, j], j = 1 .. n)]; add(mul(k), `in`(k, combinat:-choose(L, m))); return % end proc

A := Matrix(5, 5, shape = diagonal, symbol = a)

Matrix(%id = 36893491019324451524)

(1)

abs(A); [op(%)][1]; [op(`%%`)][2]; [op(`%%%`)][3]; whattype(%); whattype(`%%`); whattype(`%%%`)

a[1, 1]*a[2, 2]*a[3, 3]*a[4, 4]*a[5, 5]

 

a[1, 1]

 

a[2, 2]

 

a[3, 3]

 

indexed

 

indexed

 

indexed

(2)

komb(5, 5); [op(%)][1]; [op(`%%`)][2]; [op(`%%%`)][3]; whattype(%); whattype(`%%`); whattype(`%%%`)

a[1, 1]*a[2, 2]*a[3, 3]*a[4, 4]*a[5, 5]

 

a[1, 1]

 

a[2, 2]

 

a[3, 3]

 

indexed

 

indexed

 

indexed

(3)

If it looks like an 'índexed', if it called 'indexed', then maybe.....

simplify(abs(A)-komb(5, 5))

a[1, 1]*a[2, 2]*a[3, 3]*a[4, 4]*a[5, 5]-a[1, 1]*a[2, 2]*a[3, 3]*a[4, 4]*a[5, 5]

(4)

lprint(abs(A)); lprint(komb(5, 5))

a[1,1]*a[2,2]*a[3,3]*a[4,4]*a[5,5]
a[1,1]*a[2,2]*a[3,3]*a[4,4]*a[5,5]

 

lprint(abs(A))-lprint(komb(5, 5))

a[1,1]*a[2,2]*a[3,3]*a[4,4]*a[5,5]
a[1,1]*a[2,2]*a[3,3]*a[4,4]*a[5,5]

 

0

(5)

``

Download indexed_revisited.mw

I had this problem in the past. I guess it has something to do with an update of the Java driver. But I am not sure.

@Thomas Richard 

 

Thanks a lot. Missed that one completely!

Sorry!

1 2 3 4 5 Page 1 of 5