rkube
December 18, 2021, 7:35pm
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
Don’t use the normal equations. Just use x \ y which will use a QR factorization and be more accurate.
oheil
December 18, 2021, 7:56pm
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
2 Likes
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
