Can you prove your formula for Tensor Contractions?

[ArashMehraban]

I am writing this because I am convinced the following formulation is inaccurate mathematically. I could not verify it numerically either: The following claims: D0 = the derivative of basis functions wrt X (xi) evaluated at quadrature points D1 = the derivative of basis functions wrt Y (eta) evaluated at quadrature points D2 = the derivative of basis functions wrt Z (zeta) evaluated at quadrature points

But numerical testing shows:

D2 = the derivative of basis functions wrt X (xi) evaluated at quadrature points D1 = the derivative of basis functions wrt Y (eta) evaluated at quadrature points D0 = the derivative of basis functions wrt Z (zeta) evaluated at quadrature points

The following is a small MATLAB code that checks the above claim analytically in two ways:

1. Testing the value of D0 against the kron(B,kron(B,D))
2. Testing the value of D0 against the tensor contractions using for loops developed in HPGMG. You notice in the results of the code above that: kron(B,kron(B,D)) is equal to D0, while the formula above claims D0 = kron(D,kron(B,B))

Therefore, mathematically, I am convinced the formula above does not make sense. Now, my bigger issue is that if the tensor contractions in libCEED are implemented similar to HPGMG code (example for Poisson problem), then it appears, we are doing the exact opposite order of operations:

To be completely sure, in case I made a silly mistake in my MATLAB code, I re-implemented the tensor contractions in HPGMM and FEBasisEval from libCEED in one independent pure C code. To explicitly form kron(B,kron(B,D)) using the Tensor() function, I fed it the columns of an identity matrix and stored the results into another matrix column-wise. Please note that, the C code verifies my MATLAB results:

I have provided the code in MATLAB and C in this repo if you would like to test it yourself: https://github.com/ArashMehraban/For-libCEED-tensor-issue

To run the MATLAB code, run shape_test.m

[jeremylt]

Our test suite verifies that we can accurately interpolate and take derivatives in multiple dimensions: https://github.com/CEED/libCEED/blob/main/tests/t313-basis.c https://github.com/CEED/libCEED/blob/main/tests/t314-basis.c

[jedbrown]

@ArashMehraban This is a matter of C versus Fortran ordering conventions. When you have (I⊗A)x, does it apply A along the contiguous dimension of the data x or along the strided dimension? Matlab kron uses the Fortran convention while my paper and libCEED uses C convention.

[ArashMehraban]

@jedbrown Please note that, mathematically speaking: D @ B @ B came out to be D2 from the C code too. Would you please be able to provide me with the matrix-decomposition of your kron? (D@I@I) (I@B@I) (I@I@B)? <-- this appears to be Paul Fischer's formula but it only works for nxn matrices properly.

I have tested both numpy and Matlab the same way. They come out to be the same. Can you re-produce D0 = D @ B @ B in numpy?

[ArashMehraban]

def tensor(P,Q,name):
  4     A = np.random.rand(Q,P)
  5     B = np.random.rand(Q,P)
  6     C = np.random.rand(Q,P)
  7     Ip = np.eye(P)
  8     Iq = np.eye(Q)
  9     if name == 'Paul':
 10         IIC = np.kron(Ip,np.kron(Ip,C))
 11         IBI = np.kron(Iq,np.kron(B,Ip))
 12         AII = np.kron(Iq,np.kron(Iq,A))
 13         return np.kron(A,np.kron(B,C)) - AII.dot(IBI).dot(IIC)
 14     if name == 'Arash':
 15         IIC = np.kron(Iq,np.kron(Iq,C))
 16         IBI = np.kron(Iq,np.kron(B,Ip))
 17         AII = np.kron(A,np.kron(Ip,Ip))
 18         return np.kron(A,np.kron(B,C)) - IIC.dot(IBI).dot(AII)
 19 
 20 if __name__ == '__main__':
 21     P = 2
 22     Q = 2
 23     print('Paul for square matrices\n')
 24     print(tensor(P,Q,'Paul'))
 25     print('\nArash for nxn matrices\n')
 26     print(tensor(P,Q,'Arash'))
 27 
 28     P = 2
 29     Q = 3
 30     print('Paul for non-square matrices\n')
 31     print(tensor(P,Q,'Paul'))
 32     print('\nArash for non-square matrices\n')
 33     print(tensor(P,Q,'Arash'))

[jeremylt]

libCEED tests t313 and t314 in Python: https://replit.com/@jeremylt/FriendlyGrowingCleantech#main.py

[ArashMehraban]

@jeremylt Jeremy, I read that code in your paper already. Please note that I am not denying that the code appears to be doing the right thing. What I am asking is why can't I reproduce the mathematical formula? My approach is very simple here. I am afraid, part of the reason we are getting correct results from the tests is the fact that we use symmetric tensors.

[jeremylt]

We almost never use the same number of nodes and quadrature points though. Our tensors are not symmetric.

[ArashMehraban]

Thank you all for replying. I have always enjoyed collaborating with you and @jedbrown. However, in this particular case, I am still not convinced that the mathematical formula is correct. @jedbrown, so would you say, for the sake of argument, if I were to publish a paper that I implemented all of its code in MATLAB, should I re-introduce: D0, D1 and D2 in the oppose order of tensor products in your 2010 paper?

[jeremylt]

If you look at the repl I shared, you can see that the layout of the input needs to agree with the assumed layout for the implementation. Numpy's kron() (and MATLAB's) assume a different layout than libCEED. Our implementation is correct when you feed in correct data. I show both using np.kron() and a libCEED basis.apply() in the repl.

Whatever you write up should be consistent between how you represent mesh data, how you implement tensor contractions, and how you mathematically describe them.

Side Q: Why use MATLAB instead of libCEED? We have user friendly Python and Julia interfaces. We have performance portability. We integrate nicely with PETSc.

[ArashMehraban]

Let me answer the side Q first if I may: I use MATLAB to independently verify things that I have doubt about or want to learn their detaials. I say to myself, I try to determine if I can come up with what libCEED, Numpy, PETSc, etc. have produced. My code has come to the rescue twice already:) If you recall because of my MATLAB code we found the bug in DMPlex boundary code. Also, we were in your house and you were working on the preconditioning part of libCEED code and something was not right, then we discussed my MATLAB code and I am not sure if you recall, you said, "looking at your code, I just realized what I am multiplying extra". I was assembling a Q2 matrix but you needed its diagonal only. We ended up having scotch that day :D

Anyhow, what I am concerned about is that I cannot verify the mathematical formula with any programming language, MATLAB, Numpy, C. Don't you agree that Math must be verified anywhere? On a separate note, I also need to approach you about using libCEED and developing a better user manual. I am not complaining, I wan to make it better! I am just so busy :)

[ArashMehraban]

FWIW, I am wondering, for the sake of argument, if you changed the order of D0 and D2, if you would get more accurate results in your tests?

[jeremylt]

Ooof, that mustache on me was a bit patchy that day....

I understand why we use MATLAB to help debugging, I'm just not sure if it can be used to publish as compelling of results as libCEED can produce.

Anyhow, what I am concerned about is that I cannot verify the mathematical formula with any programming language, MATLAB, Numpy, C. Don't you agree that Math must be verified anywhere?

I agree that it should be verifiable, which is why I made the repl. In that repl, I start with f(x, y) = xx + yy + xy + x and make a vector of function values on [-1, 1]^2 . Then I used B ⊗ B, D ⊗ B, and B ⊗ D to verify that we correctly compute f, dfdx, and dfdy on quadrature points. So long as I order the input vector correctly for the tensor product/Kronecker product implementation, everything works out correctly. This is what Jed was talking about for the C vs Fortran ordering.

If you fork the repl, you can test different numbers of nodes, different numbers of quadrature points, and different functions.

On a separate note, I also need to approach you about using libCEED and developing a better user manual. I am not complaining, I wan to make it better! I am just so busy :)

I always think documentation can be better, but is there a specific area you're thinking of?

FWIW, I am wondering, for the sake of argument, if you changed the order of D0 and D2, if you would get more accurate results in your tests?

More accurate than 1E-14? My Numpy kron test has one more digit of accuracy than my libCEED basis.apply test, but I think that's because Numpy is pulling a few tricks to tighten up rounding errors because they have an assembled matrix to work with.

The libCEED integration tests have a higher error tolerance, but that's because we intentionally chose tanh because it is hard to integrate. The libCEED differentiation test has a pretty low error tolerance because it is working with a polynomial.

[ArashMehraban]

OK! That is the point! I was thinking if the formula works correctly for a certain layout, maybe we should have mentioned it! We got convergence out of nearly incompressible Neo-Hookean hyperelasticy code. So the code must be doing something right. However, I was just wondering, for the sake of mathematically sound formulas, why I cannot verify the equations. That part is not code, it is math.

As far as manuals goes, I feel like we have too complicated of examples. We need much more simple examples on how to utilize libCEED functions. It is going to be a big effort to have what I have in mind developed, but it will be worth it I think.

[jeremylt]

Have you seen the Python tutorials: https://github.com/CEED/libCEED/tree/main/examples/python

We could perhaps highlight them better in the user manual, but they have lots of smaller examples.

[ArashMehraban]

Yes! I would say, in addition to what the code does, I would elaborate on why would someone want to use such functions, in addition to coupe of use-cases. It is demonstrating what it is doing but it is not very clear why. Also, I would make so many more examples to cover more bases, especially with the C code. believe me, I constantly read the updates and examples. I have a few multigrid questions for you too. I just need to figure out a time that work for both of us :)