julialinearalgebra / genericlinearalgebra.jl Goto Github PK
View Code? Open in Web Editor NEWGeneric numerical linear algebra in Julia
Home Page: https://julialinearalgebra.github.io/GenericLinearAlgebra.jl/
License: MIT License
genericlinearalgebra.jl's People
Contributors
![dependabot[bot] avatar](https://avatars.githubusercontent.com/in/29110?v=4 "dependabot[bot]")
![femtocleaner[bot] avatar](https://avatars.githubusercontent.com/in/4123?v=4 "femtocleaner[bot]")
Stargazers
Watchers
Forkers
mfalt chrisrackauckas zgornel invenia zhangdenglinggit w111l50n eschnett ericphanson sverek purusharths juliatagbot davidanthoff sailfish009 simonbyrne standardgalactic johnfgibson wgddd andrelmfsantos playfloor dkarrasch steffenpl ranocha nbhavikatti daanhbgenericlinearalgebra.jl's Issues
Convergence issues in non-symmetric QR algorithm
I was skimming through some LAPACK related papers [1] and stumbled upon this example:
function H(n::Int)
H = zeros(2n, 2n)
for i = 1 : 2 : 2n
H[i, i+1] = 1
H[i+1, i] = 1
end
H
end
function E(n::Int)
E = zeros(2n, 2n)
for i = 1 : (n - 1)
E[2i + 1, 2i] = 1
end
E[1, 2n] = 1
E
end
my_matrix(n::Int, η::Float64 = 1e-9) = H(n) .+ η .* E(n)Maybe the choice of shifts makes this fail, even for small matrices:
> A = my_matrix(4, 1e-3);
> GenericLinearAlgebra._schur!(A, maxiter = 100000)
ERROR: ArgumentError: iteration limit 100000 reached
> A = my_matrix(4, 1e-3);
> GenericLinearAlgebra.LAPACK2.lahqr!(A)
# okNot sure I understand yet why, but might be worth figuring out what is happening.
For reference:
julia> A
8×8 Array{Float64,2}:
0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.001
1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.001 0.0 1.0 0.0 0.0 0.0 0.0
0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.001 0.0 1.0 0.0 0.0
0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.001 0.0 1.0
0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0
[1] Day, David. How the shifted QR algorithm fails to converge and how to fix it. Tech. Report 96–0913, Sandia National Laboratories, Albuquerque, NM, 1996.
fyi UUID assigned
Assigning UUID 68167eab-4858-5bc0-b60f-9fa0a04dd39f to GenericLinearAlgebra
Faulty eigenvalues from GenericLinearAlgebra when comparing to LinearAlgebra and other software
My expectation of comparing the eigenvalues produced by GenericLinearAlgebra and LinearAlgebra is that they should be the same up to some error. I do not find this to be the case for some matrices.
Consider the following session:
julia> m = [1 2 1 1 1 1; 0 1 0 1 0 1; 1 1 2 0 1 0; 0 1 0 2 1 1; 1 1 1 0 2 0; 0 1 1 1 0 2]
6×6 Array{Int64,2}:
1 2 1 1 1 1
0 1 0 1 0 1
1 1 2 0 1 0
0 1 0 2 1 1
1 1 1 0 2 0
0 1 1 1 0 2
julia> using LinearAlgebra
julia> eigvals(m)
6-element Array{Float64,1}:
0.09678807408844635
0.9999999894632876
0.9999999999999999
1.000000010536712
2.1939365664746306
4.709275359436919
julia> using GenericLinearAlgebra
julia> evals = eigvals(big.(m))
6-element Array{Complex{BigFloat},1}:
0.09678807408844671251478377594290576622468594995566734039578341791786382578280872 - 0.0im
0.526941082983575125659836351236052467530884090152227094914001065732795390546863 + 0.0im
1.000000000000000000000000000000000000000000000000000000000000000000000000000104 - 0.0im
1.473058917016424874340163648763947532469115909847772905085998934267204609453197 + 0.0im
2.193936566474630448256084556933203300255287310678896004216260729027605120603558 + 0.0im
4.709275359436922839229131667123890933520026739365436655387955853054531053613529 + 0.0im
julia> (evals[2] + evals[4])/2
1.000000000000000000000000000000000000000000000000000000000000000000000000000035 + 0.0imWe can see that two eigenvalues do not come out as being close to 1. Hence GenericLinearAlgebra appear to produce faulty eigenvalues.
This behaviour is dependent on the precision used for BigFloat as well
julia> setprecision(64)
64
julia> eigvals(big.(m))
6-element Array{Complex{BigFloat},1}:
0.0967880740884467127188 - 0.0im
0.999999999999999999458 - 4.03274504367466385843e-10im
0.999999999999999999458 + 4.03274504367466385843e-10im
1.00000000000000000011 + 0.0im
2.19393656647463044943 + 0.0im
4.7092753594369228397 + 0.0im
julia> setprecision(121)
121
julia> eigvals(big.(m))
6-element Array{Complex{BigFloat},1}:
0.096788074088446712514783775942905766097 - 0.0im
0.57659174364529533508673210001597938362 + 0.0im
0.99999999999999999999999999999999999812 + 0.0im
1.4234082563547046649132678999840206168 + 0.0im
2.1939365664746304482560845569332033081 + 0.0im
4.7092753594369228392291316671238909251 + 0.0imVarious precision numbers gives different results, 55 different precisions appear to fail to produce the correct eigenvalues in 64:256:
julia> n = 0
0
julia> for i in 64:256
setprecision(i)
if !all(0.1 .> abs.(eigvals(big.(m)) .- [0.096, 1, 1, 1, 2.19, 4.71]))
global n += 1;
end
end
julia> n
55and when there is a mismatch it is not always the same numbers as for setprecision(256) (the default), as seen in the setprecision(121) example above.
So there appears to be an issue in the the method for eigvals() introduced by GenericLinearAlgebra.
We have also checked with Mathematica that this matrix has three eigenvalues equal to 1. The above was tested on Julia 1.4.0.
This issue was brought to my attention by Severin Lüst.
Schur factorization hangs in infinite loop
Try
for i = 1 : 1000
println(LinearAlgebra.EigenGeneral.eigvals!(randn(8, 8)))
end
and soon enough it will get stuck. Seems like an issue with singleShiftQR!.
method ambiguity from mul! with LinearAlgebra in Julia 1.3
_ _ _(_)_ | Documentation: https://docs.julialang.org
(_) | (_) (_) |
_ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 1.3.0-alpha.148 (2019-08-16)
_/ |\__'_|_|_|\__'_| | Commit b5f4e87774 (0 days old master)
|__/ |
julia> using LinearAlgebra
julia> using GenericLinearAlgebra
julia> mul!(rand(3), rand(3,3)', rand(3), true, true)
ERROR: MethodError: mul!(::Array{Float64,1}, ::Adjoint{Float64,Array{Float64,2}}, ::Array{Float64,1}, ::Bool, ::Bool) is ambiguous. Candidates:
mul!(y::Union{DenseArray{T,1}, Base.ReinterpretArray{T,1,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T,1,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T,1,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}}, adjA::Adjoint{#s626,#s625} where #s625<:Union{DenseArray{T,1}, DenseArray{T,2}, Base.ReinterpretArray{T,1,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReinterpretArray{T,2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T,1,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T,1,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}, SubArray{T,2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}} where #s626, x::Union{DenseArray{T,1}, Base.ReinterpretArray{T,1,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T,1,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T,1,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}}, alpha::Union{Bool, T}, beta::Union{Bool, T}) where T<:Union{Float32, Float64} in LinearAlgebra at /home/coey/julia/usr/share/julia/stdlib/v1.3/LinearAlgebra/src/matmul.jl:118
mul!(C::Union{DenseArray{T,1}, DenseArray{T,2}, Base.ReinterpretArray{T,1,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReinterpretArray{T,2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T,1,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T,1,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}, SubArray{T,2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}} where T, adjA::Adjoint{#s24,#s19} where #s19<:(Union{DenseArray{T,2}, Base.ReinterpretArray{T,2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T,2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}} where T) where #s24<:Number, B::Union{DenseArray{T,1}, DenseArray{T,2}, Base.ReinterpretArray{T,1,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReinterpretArray{T,2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T,1,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T,1,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}, SubArray{T,2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}} where T, α::Number, β::Number) in GenericLinearAlgebra at /home/coey/.julia/packages/GenericLinearAlgebra/Lme7U/src/juliaBLAS.jl:127
Possible fix, define
mul!(::Union{DenseArray{T<:Union{Float32, Float64},1}, Base.ReinterpretArray{T<:Union{Float32, Float64},1,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T<:Union{Float32, Float64},1,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T<:Union{Float32, Float64},1,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}}, ::Adjoint{#s626,#s625} where #s625<:Union{DenseArray{T<:Union{Float32, Float64},1}, DenseArray{T<:Union{Float32, Float64},2}, Base.ReinterpretArray{T<:Union{Float32, Float64},1,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReinterpretArray{T<:Union{Float32, Float64},2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T<:Union{Float32, Float64},1,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T<:Union{Float32, Float64},2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T<:Union{Float32, Float64},1,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}, SubArray{T<:Union{Float32, Float64},2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}} where #s626, ::Union{DenseArray{T<:Union{Float32, Float64},1}, Base.ReinterpretArray{T<:Union{Float32, Float64},1,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T<:Union{Float32, Float64},1,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T<:Union{Float32, Float64},1,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}}, ::Union{Bool, T<:Union{Float32, Float64}}, ::Union{Bool, T<:Union{Float32, Float64}})
Stacktrace:
[1] top-level scope at REPL[3]:1
inverse of symmetric matrix is not guaranteed to be symmetric
Is this the right place for issues with LinearAlgebra.jl
Here is an MWE
while true
R = randn((3,3))
R[2,1] = R[1,2]; R[3,1] = R[1,3]; R[2,3]=R[3,2]
@assert R == R'
@assert inv(R) == inv(R)'
end
svd stuck in infinite loop for certain matrices
Using svd on certain matrices leads to an infinite loop:
using LinearAlgebra, GenericLinearAlgebra
m = [1 0 0 0; 0 2 1 0; 0 1 2 0; 0 0 0 -1]
svd(m) # ok
svd(BigFloat.(m)) # stuck in infinite loop
svdvals!(Float64.(m), debug=true) # also stuck, so it's not BigFloat's fault.The svd of some permutations of the matrix works:
perm = [2,3,1,4]
svd(BigFloat.(m[perm,perm])) # okFor this particular matrix, any small perturbation (e.g. m + 1e-10*rand(4,4)) seem to work fine but there are other matrices, such as
m2 = [0.4 0 0 0 0 0 0 0; 0 0.8 0 0 0 0 1 0; 0 0 -0.3 0 0 0 0 0; 0 0 0 -0.4 0 0 0 0; 0 0 0 0 0.4 0 0 0; 0 0 0 0 0 -0.3 0 0; 0 1 0 0 0 0 0.8 0; 0 0 0 0 0 0 0 -0.4]where svdvals!(m2+1e-2*rand(8,8), debug=true) often gets stuck.
The problem may be related to the code linked below. If I replace shift^2 < eps(B.dv[n1]^2) with true, the infinite loop is avoided, but I don't know if this is a proper thing to do.
GenericLinearAlgebra.jl/src/svd.jl
Lines 181 to 188 in 9f9b46f
Test of GitHub API
schur convergence with 0 eigenvalues
It seems that the _schur! algorithm gets stuck when a row and column is zero
A = [0.0 1.0 0.0;
-1.0 0.0 0.0;
0.0 0.0 0.0]
B = [0.0 0.0 0.0;
0.0 0.0 1.0;
0.0 -1.0 0.0]
eigvals(BigFloat.(A)) # Max iters reached
eigvals(BigFloat.(B)) # Max iters reachedI fixed this in an incoming PR.
This exposed another bug in _eigvals!(S::Schur{T} where the first "2x2" block in
[0.0 0.0 0.0; 0.0 0.0 1.0; 0.0 -1.0 0]resulted in 2 zero eigenvalues, which in the end returned all zero eigenvalues.
I fixed this too, however I am not familiar with the algorithm so I might have done something stupid, but it is similar to this fix: #53
Suggest to reorder arguments to rankUpdate!
I copied the style of your rankUpdate! generic in the MixedModels package then later decided to change the order of the arguments so that the Hermitian matrix is the first argument and the scalar is the last argument. In this order the argument that is modified is first, as is becoming the recommended practice, and the scalar can be given a default value in the function declaration.
Improve reductions to bi- and tridiagonal matrices
The current versions are very basic and pretty slow compared to LAPACK. They should at a minimum be blocked but we should also consider the methods described in https://books.google.dk/books?hl=en&lr=&id=bx7pAwAAQBAJ&oi=fnd&pg=PA135&dq=Sergey+V+Kuznetsov+bidiagonal&ots=U0trQWuG9j&sig=cGZDU36MrBCKNqRp4V4b1w_55M0&redir_esc=y#v=onepage&q=Sergey%20V%20Kuznetsov%20bidiagonal&f=false and reference [5]. We might even make them run in parallel and be faster vanilla LAPACK.
The right place to start testing these would be the QR factorization. The idea is the same but the QR is much simpler.
ERROR: type Schur has no field Z
(v1.1) pkg> add GenericLinearAlgebra#master
Updating registry at `C:\Users\jas\.julia\registries\General`
Updating git-repo `https://github.com/JuliaRegistries/General.git`
Updating git-repo `https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl.git`
Resolving package versions...
[14197337] ↑ GenericLinearAlgebra v0.1.0 ⇒ v0.1.0+ #master (https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl.git)
julia> using GenericLinearAlgebra
[ Info: Recompiling stale cache file C:\Users\jas\.julia\compiled\v1.1\GenericLinearAlgebra\Tm5A3.ji for GenericLinearAlgebra [14197337-ba66-59df-a3e3-ca00e7dcff7a]
julia> m = Complex{BigFloat}.(reshape(rand(Float32, 5*5),5,5))
5×5 Array{Complex{BigFloat},2}:
0.132617+0.0im 0.742358+0.0im 0.183927+0.0im 0.632733+0.0im 0.177756+0.0im
0.125482+0.0im 0.589812+0.0im 0.178626+0.0im 0.308561+0.0im 0.856793+0.0im
0.769404+0.0im 0.827643+0.0im 0.967583+0.0im 0.829904+0.0im 0.670297+0.0im
0.115544+0.0im 0.881652+0.0im 0.300636+0.0im 0.676362+0.0im 0.77726+0.0im
0.670813+0.0im 0.857607+0.0im 0.0240812+0.0im 0.29316+0.0im 0.952149+0.0im
julia> log(m)
ERROR: type Schur has no field Z
Stacktrace:
[1] getproperty(::GenericLinearAlgebra.Schur{Complex{BigFloat},Array{Complex{BigFloat},2}}, ::Symbol) at C:\Users\jas\.julia\packages\GenericLinearAlgebra\8M4H5\src\eigenGeneral.jl:72
[2] log(::Array{Complex{BigFloat},2}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.1\LinearAlgebra\src\dense.jl:655
[3] top-level scope at none:0
Add Julia matrix exponential?
I needed a pure Julia implementation of matrix exponentiation, so I updated an old 2012 gist to Julia 1.0, see here. I am happy to contribute it to this repo if desired. Probably a better implementation can be made that is closer to the current Julia implementation. After all, it only uses LAPACK.gebal! and LAPACK.gesv! as far as I can tell.
Move to JuliaMath?
BoundsError while calling svd
The following is using GenericLinearAlgebra on master:
julia> A = rand(BigFloat, 1, 2)
1×2 Matrix{BigFloat}:
0.343118 0.899983
julia> svd(A)
ERROR: BoundsError: attempt to access 1×2 Matrix{BigFloat} at index [1:2, 1:2]
Stacktrace:
[1] throw_boundserror(A::Matrix{BigFloat}, I::Tuple{UnitRange{Int64}, UnitRange{Int64}})
@ Base .\abstractarray.jl:691
[2] checkbounds
@ .\abstractarray.jl:656 [inlined]
[3] _setindex!
@ .\multidimensional.jl:917 [inlined]
[4] setindex!
@ .\abstractarray.jl:1315 [inlined]
[5] svd!(A::Matrix{BigFloat}; tol::BigFloat, full::Bool, alg::LinearAlgebra.DivideAndConquer, debug::Bool)
@ GenericLinearAlgebra C:\Users\lkape\.julia\packages\GenericLinearAlgebra\z4KPO\src\svd.jl:570
[6] svd(A::Matrix{BigFloat}; full::Bool, alg::LinearAlgebra.DivideAndConquer)
@ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.8\LinearAlgebra\src\svd.jl:176
[7] svd(A::Matrix{BigFloat})
@ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.8\LinearAlgebra\src\svd.jl:176
[8] top-level scope
@ REPL[19]:1ERROR: type Shur has no field T
julia> using GenericLinearAlgebra
julia> m = Complex{BigFloat}.(reshape(rand(Float32, 5*5),5,5))
5×5 Array{Complex{BigFloat},2}:
0.506333+0.0im 0.534346+0.0im 0.505427+0.0im 0.636836+0.0im 0.429265+0.0im
0.883622+0.0im 0.875253+0.0im 0.676688+0.0im 0.427876+0.0im 0.894883+0.0im
0.830045+0.0im 0.493769+0.0im 0.883085+0.0im 0.398726+0.0im 0.0346483+0.0im
0.00872445+0.0im 0.995922+0.0im 0.958189+0.0im 0.757107+0.0im 0.746112+0.0im
0.322169+0.0im 0.182816+0.0im 0.491731+0.0im 0.400665+0.0im 0.043752+0.0im
julia> log(m)
ERROR: type Schur has no field T
Stacktrace:
[1] getproperty at .\Base.jl:20 [inlined]
[2] log(::Array{Complex{BigFloat},2}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\LinearAlgebra\src\dense.jl:651
[3] top-level scope at REPL[3]:1
Turning on 2FA in this org
Merge `RowEchelon.jl` into `LinearAlgebra.jl`
Testing the REST API
Numerical issues with svd
I've been running into some problems with svd when multiple singular values have the same magnitude. Consider for example the following.
import GenericLinearAlgebra
# Construct a matrix with singular values [1; 1; 1; 0]
U0, _, V0 = svd(big.(reshape(0:15,4,4)))
A = U0[:, 1:3] * V0[:, 1:3]'
U, S, V = svd(A)
norm(A - U*Diagonal(S)*V')This gives an error of 2.07 (!). With GenericSVD the same error is 5.6e-77. Casting A to Float64 gives an error of 1.9e-16.
Note that the singular values themselves seem to be computed alright.
Convergence problems of some matrices
Hi,
I'm encountering some matrices that fail to converge when computing eigenvalues. Consider the following:
julia> t = [1 -2 1 -1 -1 0; 0 1 0 1 0 1; 1 -1 2 0 -1 0; 0 1 0 2 1 1; 1 0 1 0 0 0; 0 -1 1 -1 -2 0]
6×6 Array{Int64,2}:
1 -2 1 -1 -1 0
0 1 0 1 0 1
1 -1 2 0 -1 0
0 1 0 2 1 1
1 0 1 0 0 0
0 -1 1 -1 -2 0
julia> using LinearAlgebra
julia> using GenericLinearAlgebra
julia> eigvals(t)
6-element Array{Complex{Float64},1}:
-7.294730590901527e-9 + 0.0im
7.2947305672186945e-9 + 0.0im
1.4999999999999996 - 0.8660254037844369im
1.4999999999999996 + 0.8660254037844369im
1.5000000000000018 - 0.8660254037844406im
1.5000000000000018 + 0.8660254037844406im
julia> eigvals(big.(t))
ERROR: ArgumentError: iteration limit 600 reached
Stacktrace:
[1] _schur!(::GenericLinearAlgebra.HessenbergFactorization{BigFloat,Array{BigFloat,2},GenericLinearAlgebra.Householder{BigFloat,S} where S<:(Union{DenseArray{T,1}, Base.ReinterpretArray{T,1,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T,1,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T,1,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where
A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}} where T)}; tol::BigFloat, debug::Bool, shiftmethod::Symbol, maxiter::Int64, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/user/.julia/packages/GenericLinearAlgebra/sJLhE/src/eigenGeneral.jl:103
[2] _schur! at /home/user/.julia/packages/GenericLinearAlgebra/sJLhE/src/eigenGeneral.jl:91 [inlined]
[3] #_schur!#21 at /home/user/.julia/packages/GenericLinearAlgebra/sJLhE/src/eigenGeneral.jl:183 [inlined]
[4] _schur! at /home/user/.julia/packages/GenericLinearAlgebra/sJLhE/src/eigenGeneral.jl:183 [inlined]
[5] _eigvals!(::Array{BigFloat,2}; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/user/.julia/packages/GenericLinearAlgebra/sJLhE/src/eigenGeneral.jl:260
[6] _eigvals! at /home/user/.julia/packages/GenericLinearAlgebra/sJLhE/src/eigenGeneral.jl:260 [inlined]
[7] eigvals!(::Array{BigFloat,2}; sortby::typeof(LinearAlgebra.eigsortby), kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/user/.julia/packages/GenericLinearAlgebra/sJLhE/src/eigenGeneral.jl:274
[8] eigvals! at /home/user/.julia/packages/GenericLinearAlgebra/sJLhE/src/eigenGeneral.jl:271 [inlined]
[9] #eigvals#81 at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/eigen.jl:326 [inlined]
[10] eigvals(::Array{BigInt,2}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/eigen.jl:326
[11] top-level scope at REPL[33]:1
julia> GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(BigFloat.(t)); maxiter=10000))
ERROR: ArgumentError: iteration limit 10000 reached
[stacktrace]
julia> GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(BigFloat.(t)); debug=true))
Split! Subdiagonal element is: 3.144e-77 and istart now 4
block start is: 4, block end is: 6, d: 8.096e-02, t: 1.663e+00
Francis double shift! Subdiagonal is: 1.144e+00, last subdiagonal is: 8.664e-01
Split! Subdiagonal element is: 0.000e+00 and istart now 4
block start is: 4, block end is: 6, d: 1.835e-02, t: 1.380e+00
Francis double shift! Subdiagonal is: -7.721e-02, last subdiagonal is: 7.052e-01
Split! Subdiagonal element is: 0.000e+00 and istart now 4
block start is: 4, block end is: 6, d: 5.610e-06, t: 1.509e+00
Francis double shift! Subdiagonal is: -7.490e-04, last subdiagonal is: 1.118e+00
Split! Subdiagonal element is: 0.000e+00 and istart now 4
block start is: 4, block end is: 6, d: -7.790e-10, t: 1.665e+00
Francis double shift! Subdiagonal is: -3.435e-09, last subdiagonal is: 7.440e-01
Split! Subdiagonal element is: 0.000e+00 and istart now 4
block start is: 4, block end is: 6, d: -4.071e-19, t: 1.386e+00
Francis double shift! Subdiagonal is: -1.804e-18, last subdiagonal is: 7.020e-01
Split! Subdiagonal element is: 0.000e+00 and istart now 4
block start is: 4, block end is: 6, d: -9.812e-39, t: 1.469e+00
Francis double shift! Subdiagonal is: 3.850e-37, last subdiagonal is: 1.116e+00
Split! Subdiagonal element is: 0.000e+00 and istart now 4
block start is: 4, block end is: 6, d: 7.081e-76, t: 1.671e+00
Francis double shift! Subdiagonal is: -3.069e-75, last subdiagonal is: 7.506e-01
Split! Subdiagonal element is: 1.628e-150 and istart now 6
Bottom deflation! Block size is one. New iend is 5
Split! Subdiagonal element is: 0.000e+00 and istart now 4
Bottom deflation! Block size is two. New iend is 3
block start is: 1, block end is: 3, d: 1.000e+00, t: 2.000e+00
Wilkinson-like shift! Subdiagonal is: -1.000e+00, last subdiagonal is: -1.414e+00
block start is: 1, block end is: 3, d: 1.000e+00, t: 2.000e+00
Francis double shift! Subdiagonal is: -7.071e-01, last subdiagonal is: -1.000e+00
block start is: 1, block end is: 3, d: 1.000e+00, t: 2.000e+00
Francis double shift! Subdiagonal is: -1.000e+00, last subdiagonal is: 1.414e+00
block start is: 1, block end is: 3, d: 1.000e+00, t: 2.000e+00
Francis double shift! Subdiagonal is: 1.414e+00, last subdiagonal is: -7.071e-01
block start is: 1, block end is: 3, d: 1.000e+00, t: 2.000e+00
Francis double shift! Subdiagonal is: -7.071e-01, last subdiagonal is: -1.000e+00
block start is: 1, block end is: 3, d: 1.000e+00, t: 2.000e+00
Francis double shift! Subdiagonal is: -1.000e+00, last subdiagonal is: 1.414e+00
[more of the same type of debug messages]
Francis double shift! Subdiagonal is: -1.414e+00, last subdiagonal is: -7.071e-01
block start is: 1, block end is: 3, d: 1.000e+00, t: 2.000e+00
Francis double shift! Subdiagonal is: -7.071e-01, last subdiagonal is: 1.000e+00
block start is: 1, block end is: 3, d: 1.000e+00, t: 2.000e+00
Francis double shift! Subdiagonal is: 1.000e+00, last subdiagonal is: -1.414e+00
block start is: 1, block end is: 3, d: 1.000e+00, t: 2.000e+00
Wilkinson-like shift! Subdiagonal is: -1.414e+00, last subdiagonal is: -7.071e-01
ERROR: ArgumentError: iteration limit 600 reached
[stacktrace]The above seems to show that indeed it is stuck in some infinite loop.
The analytical eigenvalues are [0,0,1/2*(3-√3im),1/2*(3+√3im),1/2*(3-√3im),1/2*(3+√3im)], and this only seems to happen with non-diagonalizable matrices.
These are other matrices I have so far found that displays a similar behaviour:
[1 -2 1 -1 -1 0; 0 1 0 1 0 1; 1 -1 2 0 -1 0; 0 1 0 2 1 1; 1 0 1 0 0 0; 0 -1 1 -1 -2 0]
[1 0 -1 0 0 0; 0 1 1 1 -1 0; 1 0 -1 0 0 0; 1 -1 0 -1 1 1; 0 0 1 0 0 0; -1 0 1 0 0 0]
[1 -2 -1 1 -1 0; 0 1 0 -1 0 1; -1 1 2 0 1 0; 0 -1 -2 2 -1 -1; 1 0 -1 0 0 0; 0 -1 -1 1 0 0]
[2 0 -1 0 0 1; 0 2 0 -1 -1 0; -1 0 1 0 0 -1; 0 -1 0 1 1 0; 0 1 0 -1 0 0; -1 0 1 0 0 0]
[1 0 1 0 -1 0; -2 1 2 -1 1 1; -1 0 -1 0 1 0; 2 1 2 -1 -2 1; 1 0 1 0 -1 0; 1 -1 -2 1 0 -1]
Generic LDL' factorization
From: JuliaLang/julia#10953
@mschauer wrote:
It would be nice to have a the semidefinite Cholesky or LDL' decomposition, which allows semidefinite matrices, is more stable and avoids computation of square roots. This makes most sense as generic julia algorithm, because for BLAS floats pivoted Cholesky does the trick. I remember that pivoting was difficult for abstract matrices.
Hermitian{<:Complex} eigensolver triggers ambiguity
julia> eigen(Hermitian(big.(complex.(randn(4,4), randn(4,4)))))
ERROR: MethodError: LinearAlgebra.eigen(::Hermitian{Complex{BigFloat},Array{Complex{BigFloat},2}}) is ambiguous. Candidates:
eigen(A::Hermitian) in GenericLinearAlgebra at /Users/andreasnoack/.julia/packages/GenericLinearAlgebra/s9cSL/src/eigenSelfAdjoint.jl:630
eigen(A::Union{Hermitian{T,S}, Hermitian{Complex{T},S}, Symmetric{T,S}} where S where T<:Real) in LinearAlgebra at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/symmetric.jl:490
Possible fix, define
eigen(::Union{Hermitian{T,S} where S<:(AbstractArray{#s549,2} where #s549<:T) where T<:Real, Hermitian{Complex{T},S} where S<:(AbstractArray{#s549,2} where #s549<:Complex{T}) where T<:Real})
Stacktrace:
[1] top-level scope at none:0Add GenericSVD.jl here?
It seems like it fits nicely here, and there is no need to have it as a separate package. @simonbyrne
Need to include `using LinearAlgebra` for `svd` to work
using GenericLinearAlgebra
svd(big.([1,2;3,4]))
returns
UndefVarError: svd not defined
Stacktrace:
[1] top-level scope
@ In[2]:1
[2] eval
@ .\boot.jl:373 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base .\loading.jl:1196
If I add using LinearAlgebra before using GenericLinearAlgebra, it works.
This is, to me, a bit counter-intuitive.
SVD U is not always orthonormal
using GenericLinearAlgebra
A = zeros(BigFloat, 2, 2)
F = GenericLinearAlgebra.svd(A, full = false)
LinearAlgebra.SVD{BigFloat, BigFloat, Matrix{BigFloat}}
U factor:
2×2 Matrix{BigFloat}:
0.0 0.0
0.0 0.0
singular values:
2-element Vector{BigFloat}:
0.0
0.0
Vt factor:
2×2 Matrix{BigFloat}:
1.0 0.0
0.0 1.0U is not orthonormal, though Vt is. It would be more convenient and consistent with e.g. LAPACK if U and Vt are both always orthonormal.
juliaBLAS needs update (and tests)
The juliaBLAS.jl source file was last updated many years ago and LinearAlgebra.jl has been reorganized since. There is some inconsistency in juliaBLAS.jl with names that are now in the LinearAlgebra.BLAS module. For example, the qualified names BLAS.ger! and BLAS.syr! occur in lines 12 and 31 but unqualified names syrk! and herk! in lines 52 and 55
eigen! hangs in infinite loop
using LinearAlgebra
using GenericLinearAlgebra
M = zeros(3,3)
M[1,1] = 0.01 # Case A
# M[3,3] = 0.01 # Case B
EF = GenericLinearAlgebra.eigen!(Hermitian(M), debug=true)Cause A indefinitely executes:
QL, blockstart: 1, blockend: 3, e[blockstart]: 0.000000e+00, e[blockend-1]:-0.000000e+00, μ: 0.100000Cause B indefinitely executes:
QL, blockstart: 1, blockend: 3, e[blockstart]: NaN, e[blockend-1]:NaN, μ: NaN
I suspect this might be an easy fix.
Make test coverage decent
Find a way to use LAPACK tests
Make tridiagonal QR algorithm more precise
For some reason, my implementation requires a much smaller tolerance than LAPACK to get the same precision in the eigenvectors. I'm not sure why. Things to check
- Will a QR/QL strategy be more precise? Right now I'm only doing QR.
- Is the shift calculation avoiding canceling operations?
eigenvectors calculation in eigenSelfAdjoint
Currently, symtriUpper! is not defined and symtriLower! does not return the transformation matrix such that A = QSQ', where A is the original matrix and S is symmetric tridiagonal.
Therefore, functions such as eigQL! are initialized with a 0 x n matrix, and it skips the eigenvector computation step rather than returning them.
Implementing the simple algorithm described on the wiki Householder transformation page shows that the eigenvector calculation works as intended.
On a small 3x3 test matrix it has 65 allocations (4.063 KB) vs 17 allocations (864 bytes) for symtriLower!
I imagine adjusting the latter code would be better than adding what I've written.
EDIT:
What I wrote here was probably dumb. I strongly suspect I had just misread things.
My apologies.
Checklist ticking off implemented algorithms with references?
This is a great package! It would be nice to have a checklist with references to help contributors or generally curious people much like this issue JuliaLinearAlgebra/IterativeSolvers.jl#1.
Add Julia matrix exponential?\
svdvals! failing for empty matrices
The singular value decompositon fails for empty matrices with element types not <: BlasFloat:
julia> svdvals!(Float16.(ones(10,0)))
ERROR: DimensionMismatch("length of diagonal vector is 0, length of off-diagonal vector is 0")
Stacktrace:
[1] Bidiagonal
@ /home/julia/julia/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/bidiag.jl:11 [inlined]
[2] Bidiagonal
@ /home/julia/julia/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/bidiag.jl:17 [inlined]
[3] Bidiagonal
@ /home/julia/julia/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/bidiag.jl:75 [inlined]
[4] bidiagonalize!(A::Matrix{Float16})
@ GenericLinearAlgebra ~/.julia/packages/GenericLinearAlgebra/JoKjw/src/svd.jl:269
[5] svdvals!(A::Matrix{Float16}; tol::Float16, debug::Bool)
@ GenericLinearAlgebra ~/.julia/packages/GenericLinearAlgebra/JoKjw/src/svd.jl:487
[6] svdvals!(A::Matrix{Float16})
@ GenericLinearAlgebra ~/.julia/packages/GenericLinearAlgebra/JoKjw/src/svd.jl:487
[7] top-level scope
@ REPL[121]:1
julia> svdvals!(BigFloat.(ones(10,0)))
ERROR: DimensionMismatch("length of diagonal vector is 0, length of off-diagonal vector is 0")
Stacktrace:
[1] Bidiagonal
@ /home/julia/julia/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/bidiag.jl:11 [inlined]
[2] Bidiagonal
@ /home/julia/julia/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/bidiag.jl:17 [inlined]
[3] Bidiagonal
@ /home/julia/julia/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/bidiag.jl:75 [inlined]
[4] bidiagonalize!(A::Matrix{BigFloat})
@ GenericLinearAlgebra ~/.julia/packages/GenericLinearAlgebra/JoKjw/src/svd.jl:269
[5] svdvals!(A::Matrix{BigFloat}; tol::BigFloat, debug::Bool)
@ GenericLinearAlgebra ~/.julia/packages/GenericLinearAlgebra/JoKjw/src/svd.jl:487
[6] svdvals!(A::Matrix{BigFloat})
@ GenericLinearAlgebra ~/.julia/packages/GenericLinearAlgebra/JoKjw/src/svd.jl:487
[7] top-level scope
@ REPL[122]:1
julia> svdvals!(Float32.(ones(10,0)))
Float32[]eigen() broken for non-Hermitian matrices
Both in 1.6.2 and 1.7b3, see examples below with BigFloat.
julia> versioninfo()
Julia Version 1.6.2
Commit 1b93d53fc4 (2021-07-14 15:36 UTC)
Platform Info:
OS: macOS (x86_64-apple-darwin18.7.0)
CPU: Intel(R) Core(TM) i7-8569U CPU @ 2.80GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
julia> using LinearAlgebra,GenericLinearAlgebra
julia> eigen(Hermitian(rand(BigFloat,5,5)))
Eigen{BigFloat, BigFloat, Matrix{BigFloat}, Vector{BigFloat}}
values:
5-element Vector{BigFloat}:
-0.4887698284034353629278662405847996486390203845256632669531308786790234904305252
-0.23957178126044890148706578948820548968416904833626474092213791610934499449746
0.2420315098582791585836827261159398231227201754519835421916618918814673416784328
1.027553335067571884556797097988170714414183181354401756940097138295884750065342
2.4179948357615710761208524931932816638859780417157267213822396311464899653017
vectors:
5×5 Matrix{BigFloat}:
0.643244 0.056325 -0.703372 0.15764 0.251957
0.137664 -0.823489 0.111178 -0.380953 0.381352
0.335786 0.350451 0.593617 0.159468 0.621788
-0.395071 0.390543 -0.309542 -0.627235 0.449617
-0.546308 -0.208215 -0.211448 0.641226 0.449792
julia> eigen(rand(BigFloat,5,5))
ERROR: MethodError: no method matching eigen!(::Matrix{BigFloat}; permute=true, scale=true, sortby=LinearAlgebra.eigsortby)
Closest candidates are:
eigen!(::StridedMatrix{T}; permute, scale, sortby) where T<:Union{Float32, Float64} at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/eigen.jl:148
eigen!(::StridedMatrix{T}; permute, scale, sortby) where T<:Union{ComplexF32, ComplexF64} at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/eigen.jl:171
eigen!(::StridedMatrix{T}, ::StridedMatrix{T}; sortby) where T<:Union{Float32, Float64} at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/eigen.jl:427 got unsupported keyword arguments "permute", "scale"
...
Stacktrace:
[1] eigen(A::Matrix{BigFloat}; permute::Bool, scale::Bool, sortby::typeof(LinearAlgebra.eigsortby))
@ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/eigen.jl:237
[2] eigen(A::Matrix{BigFloat})
@ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/eigen.jl:235
[3] top-level scope
@ REPL[6]:1
julia> versioninfo()
Julia Version 1.7.0-beta3.0
Commit e76c9dad42 (2021-07-07 08:12 UTC)
Platform Info:
OS: macOS (x86_64-apple-darwin18.7.0)
CPU: Intel(R) Core(TM) i7-8569U CPU @ 2.80GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-12.0.0 (ORCJIT, skylake)
julia> using LinearAlgebra,GenericLinearAlgebra
julia> eigen(Hermitian(rand(BigFloat,5,5)))
Eigen{BigFloat, BigFloat, Matrix{BigFloat}, Vector{BigFloat}, Vector{BigFloat}}
values:
5-element Vector{BigFloat}:
-0.3013651960644129641532610658746961406964259498102245867410423605298043880153657
0.120203860860071466734019164981546073833764155606255715475030075776944488735232
0.6240177079049022904005106766910157195384973670823691563637390563798449544924699
0.7367367930040014998095951360312239998683781592003136115029612450139105705095599
2.332417706482406886949287685629486433665965295696010610051049348107239712251295
vectors:
5×5 Matrix{BigFloat}:
-0.0392047 -0.197348 0.861579 0.316458 0.34213
-0.0925792 0.961961 0.137852 -0.0203742 0.215967
-0.716137 -0.16382 -0.297161 -0.0398325 0.608621
0.360393 -0.0939928 0.133732 -0.821625 0.410281
0.5892 -0.00360311 -0.363992 0.472001 0.545487
julia> eigen(rand(BigFloat,5,5))
ERROR: MethodError: no method matching eigen!(::Matrix{BigFloat}; permute=true, scale=true, sortby=LinearAlgebra.eigsortby, jvl=false, jvr=true, jce=false, jcv=false)
Closest candidates are:
eigen!(::StridedMatrix{T}; permute, scale, sortby, jvl, jvr, jce, jcv) where T<:Union{Float32, Float64} at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/eigen.jl:163
eigen!(::StridedMatrix{T}; permute, scale, sortby, jvl, jvr, jce, jcv) where T<:Union{ComplexF32, ComplexF64} at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/eigen.jl:206
eigen!(::StridedMatrix{T}, ::StridedMatrix{T}; sortby) where T<:Union{Float32, Float64} at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/eigen.jl:478 got unsupported keyword arguments "permute", "scale", "jvl", "jvr", "jce", "jcv"
...
Stacktrace:
[1] eigen(A::Matrix{BigFloat}; permute::Bool, scale::Bool, sortby::typeof(LinearAlgebra.eigsortby), jvl::Bool, jvr::Bool, jce::Bool, jcv::Bool)
@ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/eigen.jl:276
[2] eigen(A::Matrix{BigFloat})
@ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/eigen.jl:274
[3] top-level scope
@ REPL[6]:1
TagBot trigger issue
This issue is used to trigger TagBot; feel free to unsubscribe.
If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.
If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!
Sorting order of eigvals
In julia v1.2, the order of eigenvalues are sorted differently
From NEWS.md
Eigenvalues λ of general matrices are now sorted lexicographically by (Re λ, Im λ) (#21598).
It would be great if GenericLinearAlgebra followed the same convention.
Do you intend to tag this?
this (or thereabouts) gets eigvals(::Matrix{DoubleFloat{Float64}) to work
using LinearAlgebra
import LinearAlgebra: rmul!, lmul!
include("juliaBLAS.jl")
include("householder.jl")
include("eigenGeneral.jl")Should I copy those files into a GenericLinearAlgebra subdir and credit you and this repo?
Write README
Methods overwritten
I get a lot of warnings about GenericLinearAlgebra overwriting some of its own methods, causing precompilation to be broken.
WARNING: Method definition reflectorApply!(Union{DenseArray{T, 2}, Base.ReinterpretArray{T, 2, S, A} where S where A<:Union{Base.SubArray{T, N, A, I, true} where I<:Union{Tuple{Vararg{Real, N} where N}, Tuple{Base.AbstractUnitRange{T} where T, Vararg{Any, N} where N}} where A<:(DenseArray{T, N} where N where T) where N where T, DenseArray{T, N} where N where T}, Base.ReshapedArray{T, 2, A, MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64}, N} where N} where A<:Union{Base.ReinterpretArray{T, N, S, A} where S where A<:Union{Base.SubArray{T, N, A, I, true} where I<:Union{Tuple{Vararg{Real, N} where N}, Tuple{Base.AbstractUnitRange{T} where T, Vararg{Any, N} where N}} where A<:(DenseArray{T, N} where N where T) where N where T, DenseArray{T, N} where N where T} where N where T, Base.SubArray{T, N, A, I, true} where I<:Union{Tuple{Vararg{Real, N} where N}, Tuple{Base.AbstractUnitRange{T} where T, Vararg{Any, N} where N}} where A<:(DenseArray{T, N} where N where T) where N where T, DenseArray{T, N} where N where T}, Base.SubArray{T, 2, A, I, L} where L where I<:Tuple{Vararg{Union{Int64, Base.AbstractRange{Int64}, Base.AbstractCartesianIndex{N} where N}, N} where N} where A<:Union{Base.ReinterpretArray{T, N, S, A} where S where A<:Union{Base.SubArray{T, N, A, I, true} where I<:Union{Tuple{Vararg{Real, N} where N}, Tuple{Base.AbstractUnitRange{T} where T, Vararg{Any, N} where N}} where A<:(DenseArray{T, N} where N where T) where N where T, DenseArray{T, N} where N where T} where N where T, Base.ReshapedArray{T, N, A, MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64}, N} where N} where A<:Union{Base.ReinterpretArray{T, N, S, A} where S where A<:Union{Base.SubArray{T, N, A, I, true} where I<:Union{Tuple{Vararg{Real, N} where N}, Tuple{Base.AbstractUnitRange{T} where T, Vararg{Any, N} where N}} where A<:(DenseArray{T, N} where N where T) where N where T, DenseArray{T, N} where N where T} where N where T, Base.SubArray{T, N, A, I, true} where I<:Union{Tuple{Vararg{Real, N} where N}, Tuple{Base.AbstractUnitRange{T} where T, Vararg{Any, N} where N}} where A<:(DenseArray{T, N} where N where T) where N where T, DenseArray{T, N} where N where T} where N where T, DenseArray{T, N} where N where T}} where T, AbstractArray{T, 1} where T, Number) in module GenericSVD at /home/fredrikb/.julia/packages/GenericSVD/w8oZE/src/utils.jl:5 overwritten in module GenericLinearAlgebra at /home/fredrikb/.julia/packages/GenericLinearAlgebra/l4g5a/src/qr.jl:18.
** incremental compilation may be fatally broken for this module **
WARNING: Method definition lmul!(LinearAlgebra.Givens{T} where T, Nothing) in module GenericSVD at /home/fredrikb/.julia/packages/GenericSVD/w8oZE/src/utils.jl:28 overwritten in module GenericLinearAlgebra at /home/fredrikb/.julia/packages/GenericLinearAlgebra/l4g5a/src/svd.jl:5.
** incremental compilation may be fatally broken for this module **
WARNING: Method definition rmul!(Nothing, LinearAlgebra.Givens{T} where T) in module GenericSVD at /home/fredrikb/.julia/packages/GenericSVD/w8oZE/src/utils.jl:34 overwritten in module GenericLinearAlgebra at /home/fredrikb/.julia/packages/GenericLinearAlgebra/l4g5a/src/svd.jl:6.
** incremental compilation may be fatally broken for this module **
`eigmin(::Hermitian)` doesn't get dispached to GenericLinearAlgebra.jl
julia> x = rand(BigFloat,4,4)
4×4 Array{BigFloat,2}:
0.632329 0.143291 0.140398 0.932326
0.682252 0.833668 0.170542 0.776611
0.213518 0.324609 0.610383 0.130877
0.277469 0.558703 0.215871 0.299646
julia> eigmin(Hermitian(x + x'))
ERROR: MethodError: no method matching eigvals!(::Hermitian{BigFloat,Array{BigFloat,2}}, ::UnitRange{Int64})
Closest candidates are:
eigvals!(::SymTridiagonal{#s664,V} where V<:AbstractArray{#s664,1} where #s664<:Union{Float32, Float64}, ::UnitRange) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/tridiag.jl:248
eigvals!(::Union{Hermitian{#s664,#s663}, Hermitian{Complex{#s664},#s663}, Symmetric{#s664,#s663}} where #s663<:(Union{DenseArray{T,2}, Base.ReinterpretArray{T,2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T,2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}} where T) where #s664<:Union{Float32, Float64}, ::UnitRange) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/symmetric.jl:732
eigvals!(::Hermitian; tol, debug) at /home/eric/.julia/packages/GenericLinearAlgebra/sJLhE/src/eigenSelfAdjoint.jl:566
Stacktrace:
[1] eigvals(::Hermitian{BigFloat,Array{BigFloat,2}}, ::UnitRange{Int64}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/symmetric.jl:763
[2] eigmin(::Hermitian{BigFloat,Array{BigFloat,2}}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/symmetric.jl:806
[3] top-level scope at REPL[108]:1
[4] eval(::Module, ::Any) at ./boot.jl:331
[5] eval_user_input(::Any, ::REPL.REPLBackend) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
[6] run_backend(::REPL.REPLBackend) at /home/eric/.julia/packages/Revise/XFtoQ/src/Revise.jl:1162
[7] top-level scope at none:0
[8] eval(::Module, ::Any) at ./boot.jl:331
[9] eval_user_input(::Any, ::REPL.REPLBackend) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
[10] run_backend(::REPL.REPLBackend) at /home/eric/.julia/packages/Revise/XFtoQ/src/Revise.jl:1162
[11] top-level scope at none:0Some piracy fixes it in the BigFloat case:
LinearAlgebra.eigmin(A::LinearAlgebra.RealHermSymComplexHerm{BigFloat,<:StridedMatrix}) = minimum(eigvals(A))But maybe a generic solution is desired?
sqrt of a BigFloat matrix
Any chance of getting sqrt to work?
It seems like LinearAlgebra makes assumptions about what schur returns in this function. Example:
julia> using GenericLinearAlgebra
s
julia> A = rand(BigFloat, 2, 2);
julia> sqrt(A)
ERROR: type Schur has no field values
Stacktrace:
[1] getproperty(F::GenericLinearAlgebra.Schur{BigFloat, Matrix{BigFloat}}, s::Symbol)
@ GenericLinearAlgebra ~/.julia/packages/GenericLinearAlgebra/tfuHk/src/eigenGeneral.jl:72
[2] sqrt(A::Matrix{BigFloat})
@ LinearAlgebra /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/dense.jl:865
[3] top-level scope
@ REPL[3]:1(My use case is that taking the square root of a Gram matrix is a neat way to orthogonalize a basis for a function space.)
Eigfact for SymmetricRFP
Fix debug printing in complex eigensolver
Currently, the formatting in the debug printing assumes that the values are floats so complex value causes an exception.
Add eigen
julia> A = big.(randn(10,10));
julia> eigen(A)
ERROR: MethodError: no method matching eigen!(::Array{BigFloat,2}; permute=true, scale=true, sortby=LinearAlgebra.eigsortby)
Closest candidates are:
eigen!(::SymTridiagonal{#s617,V} where V<:AbstractArray{#s617,1} where #s617<:Union{Float32, Float64}) at /Users/sheehanolver/Projects/julia-1.2/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/tridiag.jl:200 got unsupported keyword arguments "permute", "scale", "sortby"
eigen!(::SymTridiagonal{#s617,V} where V<:AbstractArray{#s617,1} where #s617<:Union{Float32, Float64}, ::UnitRange) at /Users/sheehanolver/Projects/julia-1.2/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/tridiag.jl:203 got unsupported keyword arguments "permute", "scale", "sortby"
eigen!(::SymTridiagonal{#s617,V} where V<:AbstractArray{#s617,1} where #s617<:Union{Float32, Float64}, ::Real, ::Real) at /Users/sheehanolver/Projects/julia-1.2/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/tridiag.jl:208 got unsupported keyword arguments "permute", "scale", "sortby"
...
Stacktrace:
[1] #eigen#58(::Bool, ::Bool, ::typeof(LinearAlgebra.eigsortby), ::typeof(eigen), ::Array{BigFloat,2}) at /Users/sheehanolver/Projects/julia-1.2/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/eigen.jl:139
[2] eigen(::Array{BigFloat,2}) at /Users/sheehanolver/Projects/julia-1.2/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/eigen.jl:137
[3] top-level scope at REPL[210]:1`HessenbergFactorization` does not have property `Q`
Which the LinearAlgebra version has.
julia> H = hessenberg(randn(3,3))
Hessenberg{Float64, UpperHessenberg{Float64, Matrix{Float64}}, Matrix{Float64}, Vector{Float64}, Bool}
Q factor:
3×3 LinearAlgebra.HessenbergQ{Float64, Matrix{Float64}, Vector{Float64}, false}:
1.0 0.0 0.0
0.0 -0.783288 0.621658
0.0 0.621658 0.783288
H factor:
3×3 UpperHessenberg{Float64, Matrix{Float64}}:
0.891223 0.891774 -1.97933
-0.954603 0.0322819 -1.45461
⋅ 0.259812 -0.536947
julia> H.Q
3×3 LinearAlgebra.HessenbergQ{Float64, Matrix{Float64}, Vector{Float64}, false}:
1.0 0.0 0.0
0.0 -0.783288 0.621658
0.0 0.621658 0.783288Are there a Chordal sparse matrix function?
I work in Convex Optimization, and I look for any function like a Python's Chordal Matrix Package.
Are there something like that in JuliaLang?
`Rotation` field in `Schur` struct is abstractly typed
This came up in a nanosoldier run. The R field in the Schur struct is abstractly typed, as it is missing the eltype of Rotation:
Should it read Rotation{T}?
What wrong ? ERROR: EigenGeneral not defined
What wrong ?
ERROR: EigenGeneral not defined
error in last line. All packages are instled.
_
_ _ ()_ | A fresh approach to technical computing
() | () () | Documentation: http://docs.julialang.org
_ _ | | __ _ | Type "help()" for help.
| | | | | | |/ ` | |
| | || | | | (| | | Version 0.4.0-dev+1984 (2014-12-07 23:20 UTC)
/ |_'|||__'| | Commit 052f54e (24 days old master)
|__/ | x86_64-w64-mingw32
julia> using LinearLeastSquares
ERROR: LinearLeastSquares not found
in require at loading.jl:47
julia> A = big(randn(3,3));b = A*ones(3)
3-element Array{Base.MPFR.BigFloat,1}:
-2.9446960791946807933783247790415771305561065673828125e-01
3.236218918047749248945166300472919829189777374267578125e+00
-2.56092452427787808932890811774996109306812286376953125e+00
julia> A\b # by LU
3-element Array{Base.MPFR.BigFloat,1}:
1.000000000000000000000000000000000000000000000000000000000000000000000000000017e+00
1.000000000000000000000000000000000000000000000000000000000000000000000000000017e+00
9.999999999999999999999999999999999999999999999999999999999999999999999999999827e-01
julia> qrfact(A)\b
3-element Array{Base.MPFR.BigFloat,1}:
1.000000000000000000000000000000000000000000000000000000000000000000000000000017e+00
1.000000000000000000000000000000000000000000000000000000000000000000000000000035e+00
9.999999999999999999999999999999999999999999999999999999999999999999999999999741e-01
julia> cholfact(A'A)(A'b)
3-element Array{Base.MPFR.BigFloat,1}:
1.000000000000000000000000000000000000000000000000000000000000000000000000000017e+00
1.000000000000000000000000000000000000000000000000000000000000000000000000000017e+00
1e+00
julia> vr, vi, v3 = big(randn(3)), big(randn(3)), big(randn(3));
julia> V = [vr + vi_im vr - vi_im v3];
julia> V = [vr + vi_im vr - vi_im v3]
3x3 Array{Complex{Base.MPFR.BigFloat},2}:
-2.92480246431090773473471244869870133697986602783203125e-01-1.23301741773709269689440759520948631688952445983886
71875e-01im -6.248288026283737028876430485979653894901275634765625e-01+0e+00im
1.2486111443976664059363201886299066245555877685546875e+00-4.780390824486795975367670052946778014302253723144
53125e-01im -1.437786704485576105838617877452634274959564208984375e+00+0e+00im
4.5285454219578136214607866349979303777217864990234375e-02-7.6873228736745721767498196186352288350462913513183
59375e-02im -8.374534723584889928105212675291113555431365966796875e-01+0e+00im
julia> A = real(V*diagm([1 + im, 1 - im, 1])/V)
3x3 Array{Base.MPFR.BigFloat,2}:
2.576802034195741440834721530934809151596626451040076198981114681613706386162393e+00 -2.111621855823541931575
384837701543463270383359650712117854830826270375426030325e+00
-8.559677601741757504896778861974182932186714773709164352191451088403852069777578e+00 9.78052219886619411959
11689523504923366255761234952124648193339633989847339366e+00
-4.188140715194974510238856704548191419461734044485174519602972737487350197843463e-01 1.40012483629965998591
0674299516564706036794287389179472732815371667614610927464e+00
julia> LinearAlgebra.EigenGeneral.eigvals!(copy(A))
ERROR: LinearAlgebra not defined
julia> Pkg.clone("https://github.com/andreasnoack/LinearAlgebra.jl")
INFO: Cloning LinearAlgebra from https://github.com/andreasnoack/LinearAlgebra.jl
INFO: Computing changes...
julia> LinearAlgebra.EigenGeneral.eigvals!(copy(A))
ERROR: LinearAlgebra not defined
julia> A
3x3 Array{Base.MPFR.BigFloat,2}:
2.576802034195741440834721530934809151596626451040076198981114681613706386162393e+00 -2.111621855823541931575
384837701543463270383359650712117854830826270375426030325e+00
-8.559677601741757504896778861974182932186714773709164352191451088403852069777578e+00 9.78052219886619411959
11689523504923366255761234952124648193339633989847339366e+00
-4.188140715194974510238856704548191419461734044485174519602972737487350197843463e-01 1.40012483629965998591
0674299516564706036794287389179472732815371667614610927464e+00
julia> using LinearAlgebra
ERROR: ArrayViews not found
in require at loading.jl:47
in include at boot.jl:242
in include_from_node1 at loading.jl:128
in include at boot.jl:242
in include_from_node1 at loading.jl:128
in reload_path at loading.jl:152
in _require at loading.jl:67
in require at loading.jl:51
while loading C:\Users\SAMSUNG2.julia\v0.4\LinearAlgebra\src\cholesky.jl, in expression starting on line 6
while loading C:\Users\SAMSUNG2.julia\v0.4\LinearAlgebra\src\LinearAlgebra.jl, in expression starting on line 11
julia> Pkg.update()
INFO: Updating METADATA...
INFO: Updating LinearAlgebra...
INFO: Computing changes...
INFO: No packages to install, update or remove
julia> Pkg.add("ArrayVievs")
ERROR: unknown package ArrayVievs
in wait at task.jl:51
in sync_end at task.jl:311
in add at pkg/entry.jl:319
in add at pkg/entry.jl:71
in anonymous at pkg/dir.jl:28
in cd at file.jl:30
in cd at pkg/dir.jl:28
in add at pkg.jl:20
julia> Pkg.add("ArrayViews")
INFO: Cloning cache of ArrayViews from git://github.com/JuliaLang/ArrayViews.jl.git
INFO: Installing ArrayViews v0.4.8
INFO: Package database updated
julia> Pkg.update()
INFO: Updating METADATA...
INFO: Updating LinearAlgebra...
INFO: Computing changes...
INFO: No packages to install, update or remove
julia> using LinearAlgebra
julia> A
3x3 Array{Base.MPFR.BigFloat,2}:
2.576802034195741440834721530934809151596626451040076198981114681613706386162393e+00 -2.111621855823541931575
384837701543463270383359650712117854830826270375426030325e+00
-8.559677601741757504896778861974182932186714773709164352191451088403852069777578e+00 9.78052219886619411959
11689523504923366255761234952124648193339633989847339366e+00
-4.188140715194974510238856704548191419461734044485174519602972737487350197843463e-01 1.40012483629965998591
0674299516564706036794287389179472732815371667614610927464e+00
julia> LinearAlgebra.EigenGeneral.eigvals!(copy(A))
ERROR: EigenGeneral not defined
Paul
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.