Eigen for FillArrays

Author: jishnub

src/FillArrays.jl

@@ -11,7 +11,8 @@ import Base: size, getindex, setindex!, IndexStyle, checkbounds, convert,
 
 import LinearAlgebra: rank, svdvals!, tril, triu, tril!, triu!, diag, transpose, adjoint, fill!,
     dot, norm2, norm1, normInf, normMinusInf, normp, lmul!, rmul!, diagzero, AdjointAbsVec, TransposeAbsVec,
-    issymmetric, ishermitian, AdjOrTransAbsVec, checksquare, mul!, kron, AbstractTriangular
+    issymmetric, ishermitian, AdjOrTransAbsVec, checksquare, mul!, kron, AbstractTriangular,
+    eigvecs, eigvals, eigen
 
 
 import Base.Broadcast: broadcasted, DefaultArrayStyle, broadcast_shape, BroadcastStyle, Broadcasted

src/fillalgebra.jl

@@ -562,3 +562,41 @@ end
 
 triu(A::AbstractZerosMatrix, k::Integer=0) = A
 tril(A::AbstractZerosMatrix, k::Integer=0) = A
+
+# eigen
+_eigenind(λ0, n, sortby) = (isnothing(sortby) || sortby(λ0) <= sortby(zero(λ0))) ? 1 : n
+
+function eigvals(A::AbstractFillMatrix{<:Union{Real, Complex}}; sortby = nothing)
+    Base.require_one_based_indexing(A)
+    n = checksquare(A)
+    # only one non-trivial eigenvalue for a rank-1 matrix
+    λ0 = float(getindex_value(A)) * n
+    ind = _eigenind(λ0, n, sortby)
+    OneElement(λ0, ind, n)
+end
+
+function eigvecs(A::AbstractFillMatrix{<:Union{Real, Complex}}; sortby = nothing)
+    Base.require_one_based_indexing(A)
+    n = checksquare(A)
+    M = similar(A, real(float(eltype(A))))
+    n == 0 && return M
+    val = getindex_value(A)
+    ind = _eigenind(val, n, sortby)
+    # The non-trivial eigenvector is normalize(ones(n))
+    M[:, ind] .= inv(sqrt(n))
+    # eigenvectors corresponding to zero eigenvalues
+    for (i, j) in enumerate(axes(M,2)[(ind == 1) .+ (1:end-1)])
+        # The eigenvectors are v = normalize([ones(n-1); -(n-1)]), and sum(v) == 0
+        # The ordering is arbitrary,
+        # and we choose to order in terms of the number of non-zero elements
+        nrm = 1/sqrt(i*(i+1))
+        M[1:i, j] .= nrm
+        M[i+1, j] = -i * nrm
+        M[i+2:end, j] .= zero(eltype(M))
+    end
+    return M
+end
+
+function eigen(A::AbstractFillMatrix{<:Number}; sortby = nothing)
+    Eigen(eigvals(A; sortby), eigvecs(A; sortby))
+end

test/runtests.jl

@@ -2940,3 +2940,27 @@ end
     @test triu(Z, 2) === Z
     @test tril(Z, 2) === Z
 end
+
+@testset "eigen" begin
+    sortby = x -> (real(x), imag(x))
+    @testset "AbstractFill" begin
+        @testset for val in (2.0, -2, 3+2im, 4 - 5im, 2im), n in (0, 1, 4)
+            F = Fill(val, n, n)
+            M = Matrix(F)
+            @test eigvals(F; sortby) ≈ eigvals(M; sortby)
+            λ, V = eigen(F; sortby)
+            @test λ == eigvals(F; sortby)
+            @test V'V ≈ I
+            @test F * V ≈ V * Diagonal(λ)
+        end
+        @testset for MT in (Ones, Zeros), T in (Float64, Int, ComplexF64), n in (0, 1, 4)
+            F = MT{T}(n,n)
+            M = Matrix(F)
+            @test eigvals(F; sortby) ≈ eigvals(M; sortby)
+            λ, V = eigen(F; sortby)
+            @test λ == eigvals(F; sortby)
+            @test V'V ≈ I
+            @test F * V ≈ V * Diagonal(λ)
+        end
+    end
+end