Linear algebra bug?

General Usage
1

I translating some code from octave/matlab to Julia. However, the code does not work. I have isolated the problem in the enclosed octave/julia files. The octave code produces the correct result.

The problem is in the calculation of :

G = ((A’*A)^-1)*A’ # Octave
G = inv(transpose(A)*A)*transpose(A) # Julia

My system is :
uname --all
Linux optiplex 5.4.0-42-generic #46~18.04.1-Ubuntu SMP Fri Jul 10 07:21:24 UTC 2020 x86_64 x86_64 x86_64 GNU/Linux

The problems are observed on julia-1.5.0, however it is also observed in 1.4.2.

Please advise. I assume that I have misunderstood the Julia syntax.

Octave code

#!/usr/bin/octave

A = [ 30 20 10 ;
40 30 20 ;
50 40 30 ;
60 50 40 ;
70 60 50 ]

size(A)

Q = (A’*A)^-1

G = ((A’*A)^-1)*A’

A

=> G =

3.8940e-02 2.5146e-02 1.0376e-02 -6.3477e-03 -2.1118e-02
7.3242e-04 1.4648e-03 2.4414e-04 2.9297e-03 -2.1973e-03
-6.0059e-02 -3.5156e-02 -1.0742e-02 1.5625e-02 3.9062e-02

julia code

#!/opt/julia-1.5.0/bin/julia
using LinearAlgebra

A = [ 30 20 10 ;
40 30 20 ;
50 40 30 ;
60 50 40 ;
70 60 50 ]

println(size(A))

print("A = ")
display(A)
println()

Q = inv(transpose(A)*A)
print("Q = ")
display(Q)
println()

G = inv(transpose(A)*A)*transpose(A)

print("G = ")
display(G)
println()

print("A = ")
display(A)
println()

===> G = 3×5 Array{Float64,2}:
0.0078125 0.00390625 0.00195312 -0.00195312 -0.00390625
0.0615234 0.0371094 0.0136719 -0.0117188 -0.0390625
-0.0908203 -0.0527344 -0.0166016 0.0214844 0.0566406

2

the condition number of A'*A seems to be ~infinite? That would mean, both the Julia and the Octave results are wrong.
By the way, you can use the same syntax in Julia:
G = ((A'*A)^-1)*A'
Check condition number:
using LinearAlgebra
Use more precision:
B = BigFloat.(A)
cond(B'*B)

3

To complement @AndiMD’s answer: A is singular, eg

julia> svd(A).S
3-element Array{Float64,1}:
 167.01031102734288
  10.370921393501915
   1.1114435990051559e-14

so this is a badly conditioned problem.

Do you need this for linear regression? Consider some form of regularization.

2 Likes
4

As others have observed, the matrix is badly conditioned. MATLAB actually tells you this explicitely in the warning message:

>> G = ((A'*A)^-1)*A'
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  1.015773e-18.

Also, inv is really rarely needed, and is not always more accurate than the back-solve \ operator. So, you may write this elegantly in Julia like this:

julia> (A'A)\A'
3×5 Array{Float64,2}:
  0.00820313   0.0053125   0.004375  -0.00046875  -0.0053125
  0.0635937    0.039375    0.01125   -0.0090625   -0.029375
 -0.0917969   -0.0546875  -0.015625   0.0195312    0.0546875

The results are not accurate as I and others have said, though.

2 Likes
5

Many thanks. You are both right !!

When I constructed the unit test for the function I just picked a very “regular” patter that turned out to be poorly conditioned. I assumed that the octave results were correct as I have used the function for many years without problems ( with proper well-conditioned data).

You saved me a lot of time.

Again. Many thanks.

Peter

3 Likes