How to define and evaluate the so-called 'Hadamard Product' with symbolic arrays?

Questions
1

Dear Wolfram Community:

I am struggling to make the following apparent simple products with the following symbolic arrays:

V1=VectorSymbol[v1,dim1,Reals];    V2=VectorSymbol[v2,dim2,Reals];
M1=MatrixSymbol[m1,{dim1,dim2},Reals];    M2=MatrixSymbol[m2,{dim1,dim2},Reals];
M12=M1*M2;  V1M1=V1*M1;  M2V2= V2*M2;

The idea is to operate elementwise without the summations that occur with the “Dot” product or the proliferation of dimensions that I get with “TensorProduct”.
This of course can be done with “Times” in explicit arrays, but how to do this with symbolic arrays ?

After googling a bit, this product is also called “Hadamard Product”, “Entrywise product”, “element-wise product” or “Schur product”, and according to Wikipedia, the Wolfram command is just “Times”: “*”. Reference: Hadamard Product (Wikipedia)

In particular, when I do more complicated operations like:

dot= M12.V2// TensorExpand; 
AXA= M12.PseudoInverse[M12].M12// TensorExpand; 
AX= M12.Inverse[M12]// TensorExpand;

In this cases, I get an error message:

2024-08-20.png
2024-08-20.png1413×654 41.5 KB

However I am stuck with this errors and I am wondering if there is any way out with the built in Wolfram Functions, as the documentation do not say anything about it.

I would greatly appreciate any help or suggestion about how to do this. Best regards,

Alberto Silva Ariano

(Pontificia Universidad Catolica del Peru)

2

I feel like you haven’t given us the full picture of what you’re trying to do. You mentioned “elementwise”, but none of your vectors have explicit elements specified, so, what exactly are you expecting to happen here? If I just take your main question at face value…

define and evaluate the so-called ‘Hadamard Product’ with symbolic arrays

Then you’ve already figured it out. It’s just Times

M1*M2

The output isn’t particularly interesting, but without elements, I’m not sure what else you expect. As for “evaluating” it, what further evaluation do you expect to be performed?

I know you said specifically “symbolic arrays”, but is it clear to you how Times works with lists?

{a, b, c}*{d, e, f}
(* {a d, b e, c f} *)

That’s all I can think of to answer the evaluation question.

3

I want to clarify what I mean by “evaluation”. My problem is not with expressions like:

A*B

As the result is trivial, and I get no error messages at that step. The problem is with more complicated expressions like:

(A*B).C;   (A1*B1).(A2*B2);   
(M1*M2).Inverse[M1*M2];
(A1*A2).PseudoInverse[A1*A2].(A1*A2);

Where the “dot” is not evaluated and I get this error message:
Product of nonscalar expressions encountered in __

I would greatly appreciate any suggestion on how to do this second set of operations.
As for the usefulness of this type of expressions, it is very hard to compute those expressions explicitly when the arrays are very big.

CASE STUDY: chemical analysis and sampling stages

A typical chemical analysis (like an ICP assay, i.e 30-40 chemical elements per assay) of many samples (> 10 samples):

EXAMPLE: let m[i] be the mass of the i-th sample, and a[i,j] the concentrations of the j-th element in the i-th sample. The mean composition of the samples (i.e. an estimate of the lot composition) would be the solution to something like:

aLOT[j]*mLOT== mSAMPLES[i].aSAMPLES[i,j]

So far, no problem, as mLOT is a scalar.

But if the samples were taken from a set of groups that are mixed together during sampling (i.e. the sampling increments that made the samples) then the system become something like:

 mSAMPLES[i]*aSAMPLES[i,j] == nGROUPS[j,k].(mGROUPS[k]*aGROUPS[k,j])

The masses are known, and so the sample concentrations. The group concentrations are in principle unknown, but constrained (the concentrations of some groups are known or assumed) With more complicated sampling stages, the system becomes more and more complicated.

Solving it with with explicit arrays would be tough. I am hoping to solve the symbolic expression to get a simpler result where I could just replace the explicit arrays and then the data…

4

But what do you mean? Let’s just stick with one example for clarity:

(A*B).C

What do you want that to evaluate to? If I were to do that on paper, well, I’d just write

(A*B).C

To say “the dot is not evaluated” implies that you have in mind some expression that it should be evaluated to. What is that expression?

5

Let A1, B and C be known arrays. I want to solve the equation for B.
If there were no “times” inside “dot”, but just “dot”, like:

(A1.A2).B == C

This would be just:

A1.A2.B == C

And the solution would be in compact form:

A2==PseudoInverse[A1].C.PseudoInverse[B]

But how to solve this for A2? I have no idea:

(A1*A2).B == C

I am stuck here:

A1*A2==C.PseudoInverse[B]

As I said in my last (updated) comment, with more complicated systems (and big arrays) I am in trouble with explicit expressions like the ones solved by the functions Reduce[] and Solve[] …

6

Okay, I think I understand the context better. I’m wondering why you think “solving it with with explicit arrays would be tough”. I don’t know that it isn’t tough, I’m just wondering if you’ve actually tried it and determined that to be infeasible.

But as for being stuck on solving for A2 in

A1*A2 == C.PseudoInverse[B]

why can’t you just divide?

A2 == C.PseudoInverse[B] / A1

I may be way over my head here, so sorry if I’m missing something. Maybe I still don’t understand what kind of workflow you’re trying to build.

7

Are the error messages important? Are the results wrong?

The messages themselves are mysterious in that they do not arise from TensorRank:

Trace[
 AXA = M12 . PseudoInverse[M12] . M12 // TensorExpand,
 tr_TensorRank :> 
  Check["Success"[HoldForm[tr] -> tr], HoldForm[$Failed@tr]],
 TraceInternal -> True]

Every evaluation is a “Success”.

It turns out that the message comes from an internal pattern-check function used in the tensor-expand tool SymbolicTensors`SymbolicTensorsDump`$TensorExpandRules. The check prevents a rule from being applied to a product x*b if x is not a scalar. This is a good thing. I’m not sure why it deserves a message when x is not a scalar. It is possible that TensorExpand does not deal with Hadamard products. (Are there rules/identities? Not my area.) The list of rules is long, but I did not find a rule for dealing with them. I may have missed them.

If there are identities, you may be able to apply them, if TensorExpand[] does not.

8

You cannot just “divide” an array with another one. The operation A/B can de defined as the following:

A/B ==A*(B^-1)

The problem is that you cannot take powers of arrays as if they were just scalars. Imagine that some elements are equal to zero: you can get “Infinity” expressions for those elements. The Wolfram Language issues a warning: “power of nonscalar expressions is undefined”.

With arrays, you must use functions like “Inverse” and “MatrixPower” (for square matrices) and “PseudoInverse” (for general rectangular matrices). However, when computing something like this:

PseudoInverse(A*B)

The Wolfram Language sends another error message: “product of nonscalar expressions” and evaluation is stopped.

9

Thank you for your answer.

How did you got to find the origin of the error message?

You used a function called “Trace[]”, a function unknown to me, and it may be very useful to get to the core of the problem. Obviously it is not the matrix/array trace (Wolfram command for it is “Tr[]”)

10

Okay, I’m sorry for being simplistic. You’re wanting to solve things symbolically, but you only are doing this to reduce the computational burden (at least that’s what I thought you said). At the end of this whole exercise, you are presumably going to plug in actual data, as in actual numeric matrices. All I’m saying is that if you have this expression

A1 * A2 == Z

that’s just a shorthand for a whole system of equations. So, if you have numeric values for A1 and Z, then you can just solve for A2 directly. If A1 has zeros, then you can’t get solutions for those entries in your matrix anyway. If you divide, you’ll end up with Indeterminate expressions at those locations, but you’ll have values everywhere else. Again, maybe it’s just my ignorance, but I don’t understand how you can get better than that in this case. Or you could do this instead if you like:

Solve[A1 * A2 == Z]

In this formulation, you’ll need to have a full representation of A2, but that’s not hard.

Or, maybe we can just take a completely different tack. Do you know the mathematical formulas for what you want (as in how to compute the numeric matrices, not just symbolic expressions)? If so, just provide those and we can probably figure out how to turn them into WolframLanguage expressions.

Or, yet another alternative, just solve the whole system with symbolic entries to come up with a formula, and then apply that formula to your data. I would think that would be more computationally intense than just solving numerically, but at least you’d only have to do it once and all future numeric calculations would be easier (I assume).

11

Trace[] and relatives trace the steps of evaluation and are documented here: Trace—Wolfram Documentation and the “See Also” section near the bottom of the page.

TraceInternal I’ve known about for years and years. It’s not in the documentation, and I don’t know how I learned about it. Maybe an old Usenet group. Setting it to True shows more of the internal steps, which can involve megabytes or even hundred of megabytes of output. Some internal functions are not traceable (for instance, calls to the LAPACK library functions). Some function reveal some of the internal machinations with the setting False (for instance Plot). The normal usage by regular users is to debug their own programs or to understand the evaluation procedure.

How did I find my answer? Not easy to answer. First, I’ve used Mathematica for over 30 years, and I feel I can predict fairly well what types of things might be happening under the hood, so to speak. Second, I’ve used Trace a lot, and I’m familiar with the structure of the output and know some tricks. Like Control-Period for moving up the evaluation chain (or double-click and drag). Third, there are some things that are relatively new that help: GeneralUtilities`PrintDefinitions (will print the code of some of the internal functions),Position[] (to search a very large trace output), the Find menu command, and the 3-dot button on messages.

Now, in your case, the 3-dot button on the message will show the stack trace, which contains the offending line of code. Tracing evaluation for that failed, and I kept making my pattern more general and included TraceAbove to figure out how the line was called. I used Find and Control-Period, which led me to the tensor rules. Usually people who do this are befuddled the first time, and they’re surprised it takes me only a few minutes. Anyway, this gives the code in which you can find the answer, if you want to play with it:

Trace[
 AXA = M12 . PseudoInverse[M12] . M12 // TensorExpand,
 HoldPattern[If[_, _Message]],
 TraceInternal -> True, TraceAbove -> True]

But don’t get your hopes up too much. I tried to describe the core of the problem inmy previous answer. But to elaborate a little, the check keeps from being applied a rule that factors out a scalar from Dot: (x*a).b) -> x*(a.b) is valid when x is a scalar but not if it’s a tensor of rank > 0. If you pore over the tensor-expand rules, you may find that there are none that deal with Hadamard products. But I think the error messages do not indicate an error in your code.