Weird matrix multiplication results from MKL

New to Julia
, , ,
1

Hi, I’m trying to use a simple OLS solution to an overdetermined linear system using MKL and LinearAlgebra, with coefficient matrix

‘’’
M = [200000 0 0 (2e5)^2 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 2e5 0 0 (2e5)^2; 0 0 3e5 0 (3e5)^2]
‘’’
and observation

‘’’
y = [2e5 0 0 0 4e5 9e5]
‘’’
When I try to get the ols solution of the system with
‘’’
(M’ * M)^(-1) * M’ * y’
‘’’
the output is
‘’’
5×1 Matrix{Float64}:
1.9780612128670327
3.969978473845913
2.0318402149168833
9.652545605920437e-6
5.091301176166074e-6
‘’’
but this is incorrect. In particular, if I first store the result
‘’’
u = (M’ * M)^(-1) * M’
‘’’
and then compute

‘’’
u * y’
‘’’
the output is instead
‘’’
5×1 Matrix{Float64}:
1.3066744941170327
0.645027301970913
1.7964886524168833
1.436876312377322e-5
-5.193294015312655e-6
‘’’
Any pointers on what I am doing wrong/why this is happening?

2

Perhaps,

julia> det(M'M)
0.0

Might be the issue. Why is MKL allowing the inverse to be taken? (The column 4 is twice that of column 1.)

1 Like
3

Aside from the issue with the singular matrix, you almost never want to use the normal form solution (M'M)^-1 M' y, and you should use M \ y instead which uses better algorithms. The normal form doubles the condition number of the matrix to be inverted.

4 Likes