I am working on calculating pcov from a Jacobian matrix by hand, and understand most of it but am struggling to understand how the .dot function works from numpy in python. I have looked at the manual and attempted to solve it by hand, but it does not match the output from python
Here is the original Jacobian matrix:
Jacobian: [[ 1. 1. 1. ]
[ 4. 2. 0.99999998]
[ 9.00000004 3. 1.00000004]
[ 16. 4. 1.00000004]
[ 25.00000004 5.00000011 0.99999998]]
and after this line of code:
precov = jacobian.T.dot(jacobian)
this is the output:
precov = [[979.00000266 225.00000313 55.00000039]
[225.00000313 55.00000113 15.00000023]
[ 55.00000039 15.00000023 5.00000007]]
From what I understand, I am expecting the above matrix to be the dot product of the original matrix and the transposed matrix, but when I calculate this by hand I get:
[[3 7 13]
[7 21 43]
[13 43 101]]
Does anyone have any insight into what I am missing or not understanding? Thankyou!
To a close approximation, here is your original matrix.
import numpy as np
A = np.arange(1, 6)[:, np.newaxis] ** np.arange(3)[::-1]
A
# array([[ 1, 1, 1],
# [ 4, 2, 1],
# [ 9, 3, 1],
# [16, 4, 1],
# [25, 5, 1]])
When you calculate by hand, you must be doing something like:
A[:3, :3].dot(A[:3, :3].T)
# array([[ 3, 7, 13],
# [ 7, 21, 43],
# [13, 43, 91]])
That is, you seem to be taking only the first three rows of the matrix, say $B$, then computing $B B^T$.
NumPy is providing the correct result for $A^T A$.
Just look at the top left entry of the result. If you transpose your matrix, then the first row is:
A.T[0]
# array([ 1, 4, 9, 16, 25])
If you take the matrix dot product with A, then you are looking at the dot product of this vector with itself (that is, multiplying by its transpose).
By hand, that would be 1 + 16 + 81 + 256 + 625 = 979, which agrees with the NumPy result. You are getting 3, which seems to be the dot product of three ones (the first row) with itself. Note that matrix multiplication is not commutative: $A^T A != A A^T$ in general.
It may be helpful to brush up on matrix multiplication, which is what the dot method does for two 2D arrays. (For readability, consider using the @ operator instead: A.T.dot(A) == A.T @ A.)