Harry Garst

284 Reputation

6 Badges

20 years, 53 days

MaplePrimes Activity


These are replies submitted by Harry Garst

@janhardo 

That’s fantastic!
I’m doing something else right now, but I’ll dive into your work later.

Hartelijke groet,
Harry

Thanks for all the good work.

Even without being asked, this thread is moving in the right direction. What about introducing a new Boolean rule-based operator that checks whether part of an expression on the left-hand side can be replaced by a corresponding mathematical part on the right-hand side?

Below are a few rules that could serve as a starting point. Regardless of the size of the problem, the operator would return a true or false value depending on whether a successful substitution is possible, completely bypassing Maple’s usual heavy machinery.

Let’s call this new operator evilb instead of evalb. :-)
My background is in psychology, and I lack formal mathematical training. Linear algebra is quite abstract for someone like me. Also lacking a proper sense of modesty, I’m going to tell my friend that I’m now into abstract linear algebra.:-)

@dharr 

But the question you asked has more to do with subsindets, which scans through an expression and applies the inner procedure whenever it comes across "Determinant"

If there were a type called “Determinant of Kronecker Product”, or if the user were able to define such a type themselves, then a single command would be enough to solve the problem. Nice — a first step toward implementing abstract linear algebra.

@janhardo 

Thank you for your worksheet. 
I still do not get why you have to assign two new variables.
Without both newly created variables the command:
 Determinant(A)^4;
gives the same output.

Hartelijke groet,
Harry

@dharr 

I learned something new today. evalb returns false, but verify behaves as expected. It’s always a mystery to me how these Boolean operators work.

On the downside, it doesn’t give an unevaluated solution. I don’t have the courage to try dim = 4.

The procedure was proposed by ChatGPT. I had never seen a procedure inside another procedure before, and I still don’t understand how the argument d in the inner procedure call receives the return value of the outer procedure.

Thanks for your support.

Dharr_reply.mw

I could not let it go and came up with this one:

restart; with(LinearAlgebra); dim := 2; A := Matrix(dim, dim, symbol = a); B := Matrix(dim, dim, symbol = b)

expr := (''Determinant'')((''KroneckerProduct'')(''A'', ''B''))

('LinearAlgebra:-Determinant')(('LinearAlgebra:-KroneckerProduct')('A', 'B'))

(1)

if `and`(evalb(op(0, convert(indets(expr), list)[nops(convert(indets(expr), list))-1]) = Determinant), evalb(op(0, convert(indets(expr), list)[-1]) = KroneckerProduct)) then result := ('Determinant')([op([op(expr)][1])][1])^(eval(('RowDimension')([op([op(expr)][1])][1])))*('Determinant')([op([op(expr)][1])][2])^(eval(('RowDimension')([op([op(expr)][1])][2]))) end if; result; evalb(eval(expr) = expand(result))

LinearAlgebra:-Determinant(A)^2*LinearAlgebra:-Determinant(B)^2

 

(a[1, 1]*a[2, 2]-a[1, 2]*a[2, 1])^2*(b[1, 1]*b[2, 2]-b[1, 2]*b[2, 1])^2

 

true

(2)

restart; with(LinearAlgebra); dim := 25; A := Matrix(dim, dim, symbol = a); B := Matrix(dim, dim, symbol = b)

expr := (''Determinant'')((''KroneckerProduct'')(''A'', ''B''))

('LinearAlgebra:-Determinant')(('LinearAlgebra:-KroneckerProduct')('A', 'B'))

(3)

if `and`(evalb(op(0, convert(indets(expr), list)[nops(convert(indets(expr), list))-1]) = Determinant), evalb(op(0, convert(indets(expr), list)[-1]) = KroneckerProduct)) then result := ('Determinant')([op([op(expr)][1])][1])^(eval(('RowDimension')([op([op(expr)][1])][1])))*('Determinant')([op([op(expr)][1])][2])^(eval(('RowDimension')([op([op(expr)][1])][2]))) end if

LinearAlgebra:-Determinant(A)^25*LinearAlgebra:-Determinant(B)^25

(4)

Don't do this!!!

 

NULL

Obviously, you first have to add the quotes to the expression and instead you could as well write the Zehfuss determinant.

Download Attempt_256.mw

@dharr 

The problem is that Maple calculates the determinant, but in this case there is no reason for doing so. The question is whether there is an equality or not. The result of Zehfuss is general and does not depend on the sizes of the matrices.
Is there a way to detect whether there is a determinant of a Kronecker Product in an expression? If so, we can replace it.
ChatGPT offered this solution (and many others), but it did not return |A|^n * |B|^n.

ZehfussReplace := proc(expr)
    subsindets(expr, specfunc(anything, Determinant),
        proc(d)
            local arg, A, B, m, n;

            arg := op(1, d);

            if type(arg, specfunc(anything, KroneckerProduct)) then
                A := op(1, arg);
                B := op(2, arg);

                m := LinearAlgebra:-RowDimension(A);
                n := LinearAlgebra:-RowDimension(B);

                return 'Determinant'(A)^n * 'Determinant'(B)^m;
            else
                return d;
            end if;
        end proc
    );
end proc:

@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:

1 2 3 4 5 Page 1 of 5