@janhardo 

I will be happy to share my code, but first I have to prepare it so that it is understandable. This will take some time. I think I will present the code step by step; that may be the best way to explain the results to the reader in the final article I would like to publish.

Groet, Harry

@janhardo 

I was curious whether I was the only person interested in using Maple to understand factor models. I have never used LinearSolve, but I know that the parameters of confirmatory factor models can be estimated iteratively using the Newton–Raphson algorithm (or the Expectation–Maximization algorithm). Although there are some examples in the literature that demonstrate how to solve the estimation problem using non‑iterative methods, these approaches are not used in practice.

My current interest is in estimating the standard errors of the parameters of some small confirmatory factor models. This involves inverting the asymptotic covariance matrix, which—even for small models—is a huge and complicated matrix. Thanks to Maple’s powerful factorization capabilities, I was finally able to find a symbolic solution for a model that required inverting a very complicated 13×13 matrix.

Now I would like to publish the results. But how should I explain them to the reader? I want to show how the results can be derived, rather than confronting the reader with a solution that can only be obtained using Maple (and whose correctness can only be verified again by running Maple).
groet Harry

@janhardo 

I'm trying to figure out how your program works, but I suspect that you prepare the expression for use in Maple’s LinearSolve engine. The core of your program is:

X_vec := LinearAlgebra:-LinearSolve(K, vec_rhs);

Instead of using the devec operator, you use the following loop:

for i to m do
    for j to n do
        X_sol[i, j] := X_vec[(j - 1)*m + i];
    end do;
end do;

You are using some tricks I’m not yet familiar with, so it will take some time to completely unravel your code.
groet Harry

ps. somewhere I found the term factor analysis in your work. Did you work on these models?

@janhardo

Here you also use inert functions and completely bypass Maple’s linear‑algebra module. This is completely new to me and very interesting. I will study it closely.
groeten Harry

snippet.mw

@janhardo 

I have been working for more than a year on symbolically inverting the asymptotic covariance matrix for some small factor models. I know that sounds unusual—why would anyone want to do that?—but it actually makes sense if you want to understand how standard errors depend on the estimation of the parameter values. This project only makes sense if the results can be factorized into relatively simple expressions, and in my case they can. Thanks to Maple’s extraordinary capabilities for factoring large symbolic expressions, I was finally able to obtain all the results I needed.

However, in order to publish the results, I am looking for a way to guide the reader toward the final expressions. Unfortunately, splitting the problem into smaller pieces by using inverses of patterned matrices does not yield the elegant factorizations that exist for the inverse and determinant of the full matrix. Therefore, I am still looking for ways to explain the final results in an understandable way.

Perhaps your module can help by splitting the problem into smaller, more comprehensible parts. It seems you uploaded a list of features, but not the module itself.

@janhardo 

More than ten years ago, I began trying to unravel Hessian matrices in statistical models. The first stumbling blocks were how to code all the specific permutation and selection matrices, such as the duplication, commutation, and elimination matrices. I also needed code for the vec and vech operators. It took some time before I could continue my research.

It would have been very helpful if there had been a library where all those procedures were available. I know Maple is not used by many statisticians—my colleagues preferred working in R—but for understanding how the mathematics behind all this works, I sometimes prefer using symbolic matrices.

Did you know that when I first started using Maple, even the Kronecker product was missing?

groeten Harry

@janhardo 

Thanks!
Are you going to upload it to the Application section on the Maplesoft website?

groeten Harry

@janhardo 

My confusion is about this particular rule:

If matrix A is conformable with X as well as with the Kronecker product X⊗Y, then the number of rows of X must be the same as the number of rows of X⊗Y. The same reasoning applies to Y and B. Clearly, the equality only makes sense when some of the objects involved are vectors; if all of them are matrices, the expression cannot be conformable.
Do you agree?

toch.mw

@janhardo 

Today I am testing your module. I am confused about this one:

A := Matrix(2, 2, symbol = a);
B := Matrix(2, 2, symbol = b);
X := Matrix(2, 2, symbol = x);
Y := Matrix(2, 2, symbol = y);
oldexpres := (KroneckerProduct(A, B)) . (KroneckerProduct(X, Y));
newexpres := KroneckerProduct(A . X, B . Y);
simplify(oldexpres - newexpres);
Equal(oldexpres, expand~(newexpres));
                          [0  0  0  0]
                          [          ]
                          [0  0  0  0]
                          [          ]
                          [0  0  0  0]
                          [          ]
                          [0  0  0  0]

                              true

@janhardo 

en dat allemaal op Koningsdag!

@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