Strange error when performing linear regression

General Usage
1

What is going on here?

julia> x = log.(Nz_list[1:end])
5-element Vector{Float64}:
 4.1588830833596715
 4.852030263919617
 5.545177444479562
 6.238324625039508
 6.931471805599453

julia> 

julia> y = log.(dist_list[1:end])
5-element Vector{Float64}:
 -4.074961466558573
 -5.1239528564742365
 -6.165969558863518
 -7.206262977477147
 -8.246126826893256

julia> (x' * x) \ (x' * y)
-1.1233909121590568

julia> x = log.(Nz_list)
1×5 Matrix{Float64}:
 4.15888  4.85203  5.54518  6.23832  6.93147

julia> y = log.(dist_list)
1×5 Matrix{Float64}:
 -4.07496  -5.12395  -6.16597  -7.20626  -8.24613

julia> (x' * x) \ (x' * y)
ERROR: SingularException(4)
Stacktrace:
 [1] checknonsingular
   @ /builddir/build/BUILD/julia-1.6.4/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/factorization.jl:19 [inlined]
 [2] checknonsingular
   @ /builddir/build/BUILD/julia-1.6.4/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/factorization.jl:21 [inlined]
 [3] #lu!#136
   @ /builddir/build/BUILD/julia-1.6.4/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:85 [inlined]
 [4] #lu#140
   @ /builddir/build/BUILD/julia-1.6.4/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:273 [inlined]
 [5] lu (repeats 2 times)
   @ /builddir/build/BUILD/julia-1.6.4/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:272 [inlined]
 [6] \(A::Matrix{Float64}, B::Matrix{Float64})
   @ LinearAlgebra /builddir/build/BUILD/julia-1.6.4/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/generic.jl:1136
 [7] top-level scope
   @ REPL[89]:1
2

Don’t use the normal equations. Just use x \ y which will use a QR factorization and be more accurate.

3

Not an expert of Linear Algebra here, but if I remember correct:
(x' * x) is a scalar in the case of x being a vector, but a matrix in the case of x being a matrix:

julia> x=reshape([ 4.15888,  4.85203,  5.54518,  6.23832,  6.93147], 1, 5 )
1×5 Matrix{Float64}:
 4.15888  4.85203  5.54518  6.23832  6.93147

julia> x' * x
5×5 Matrix{Float64}:
 17.2963  20.179   23.0617  25.9444  28.8272
 20.179   23.5422  26.9054  30.2685  33.6317
 23.0617  26.9054  30.749   34.5926  38.4362
 25.9444  30.2685  34.5926  38.9166  43.2407
 28.8272  33.6317  38.4362  43.2407  48.0453

versus:

julia> x=[ 4.15888,  4.85203,  5.54518,  6.23832,  6.93147]
5-element Vector{Float64}:
 4.15888
 4.85203
 5.54518
 6.23832
 6.93147

julia> x' * x
158.549411991

Is this good enough for an explanation? Well, a real expert is replying… lets wait :wink:

2 Likes
4

Exactly. In the first case where x and y are vectors, x' * x and y' *y are scalar valued inner products, and (x' * x) \ (x' * y) is scalar division.

julia> x = [4.1588830833596715, 4.852030263919617, 5.545177444479562, 6.238324625039508, 6.931471805599453]
5-element Vector{Float64}:
 4.1588830833596715
 4.852030263919617
 5.545177444479562
 6.238324625039508
 6.931471805599453

julia> y = [-4.074961466558573, -5.1239528564742365, -6.165969558863518, -7.206262977477147, -8.246126826893256]
5-element Vector{Float64}:
 -4.074961466558573
 -5.1239528564742365
 -6.165969558863518
 -7.206262977477147
 -8.246126826893256

julia> x' * x
158.54949459300647

julia> x' * y
-178.113061353195

julia> (x' * x) \ (x' * y)
-1.1233909121590568

In the second case where x is a 1 x 5 matrix, x' * x is a 5x1 x 1x5 outer product matrix, which is necessarily rank-deficient. Let A = x' * x. Then the jth column of A is the vector x' times the scalar x[j]. So all the columns of A are linearly dependent, and A has rank 1.

julia> x = x'
1×5 adjoint(::Vector{Float64}) with eltype Float64:
 4.15888  4.85203  5.54518  6.23832  6.93147

julia> y = y'
1×5 adjoint(::Vector{Float64}) with eltype Float64:
 -4.07496  -5.12395  -6.16597  -7.20626  -8.24613

julia> A = x' * x
5×5 Matrix{Float64}:
 17.2963  20.179   23.0617  25.9445  28.8272
 20.179   23.5422  26.9054  30.2685  33.6317
 23.0617  26.9054  30.749   34.5926  38.4362
 25.9445  30.2685  34.5926  38.9167  43.2408
 28.8272  33.6317  38.4362  43.2408  48.0453

julia> A[:,1]
5-element Vector{Float64}:
 17.296308501055247
 20.179026584564458
 23.061744668073665
 25.944462751582876
 28.827180835092083

julia> x' * x[1]
5-element Vector{Float64}:
 17.296308501055247
 20.179026584564458
 23.061744668073665
 25.944462751582876
 28.827180835092083

julia> using LinearAlgebra

julia> rank(A)
1
1 Like