Add the Jordan canonical form

[MikaelSlevinsky]

This enables

julia> using MatrixFactorizations

julia> A = [0 0 1 7 -1; -5 -6 -6 -35 5; 1 1 -7 7 -1; 0 0 0 -9 0; 2 1 -5 -42 -3]
5×5 Matrix{Int64}:
  0   0   1    7  -1
 -5  -6  -6  -35   5
  1   1  -7    7  -1
  0   0   0   -9   0
  2   1  -5  -42  -3

julia> λ = [-9, -7, -7, -1, -1//1]
5-element Vector{Rational{Int64}}:
 -9
 -7
 -7
 -1
 -1

julia> V, J = jordan(A, λ)
Jordan{Rational{Int64}, Matrix{Rational{Int64}}, Matrix{Rational{Int64}}}
Generalized eigenvectors:
5×5 Matrix{Rational{Int64}}:
 0  0   0  1   0
 0  1   0  0  -1
 0  1  -1  0   0
 1  0   0  0   0
 7  1  -1  1  -1
Jordan normal form:
5×5 Matrix{Rational{Int64}}:
 -9   0   0   0   0
  0  -7   1   0   0
  0   0  -7   0   0
  0   0   0  -1   1
  0   0   0   0  -1

julia> A*V == V*J
true

Note that Jordan blocks are unstable with respect to perturbations in the matrix, so the code is generally best used with rational arithmetic (with matrices with rational eigenvalues), though general support is provided.

In the future, the Jordan normal form of the matrix could benefit from structuring as a BlockDiagonal matrix with a new type of LayoutMatrix to describe a JordanBlock, but this comes with the drawback of including other packages in the requirements. exp(J::JordanBlock) (and other matrix functions) is also an UpperTriangularToeplitz matrix, though this gets quite esoteric.

[codecov]

Codecov Report

❌ Patch coverage is 99.40828% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 93.79%. Comparing base (cbc6178) to head (892043b).
⚠️ Report is 5 commits behind head on master.

Files with missing lines	Patch %	Lines
src/jordan.jl	99.40%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master      #76      +/-   ##
==========================================
+ Coverage   93.17%   93.79%   +0.62%     
==========================================
  Files          10       11       +1     
  Lines        1524     1693     +169     
==========================================
+ Hits         1420     1588     +168     
- Misses        104      105       +1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[dlfivefifty]

Can't we do:

const JordanBlock{T} = Bidiagonal{T, FillVector{T}}

JordanBlock(λ, n) = Bidiagonal(Fill(λ,n), Fill(one(λ), n-1), :U)

[MikaelSlevinsky]

We could. Note also that a simple struct:

struct JordanBlock{T} <: LayoutMatrix{T}
    λ::T
    n::Int
end

@inline size(J::JordanBlock) = (J.n, J.n)
function getindex(J::JordanBlock{T}, i::Integer, j::Integer) where T
    if i == j
        return J.λ
    elseif i == j-1
        return one(T)
    else
        return zero(T)
    end
end

function rowsupport(_, J::JordanBlock, j)
    if j == J.n
        return j:j
    else
        return j:j+1
    end
end

function colsupport(_, J::JordanBlock, j)
    if j == 1
        return j:j
    else
        return j-1:j
    end
end

describes the blocks without changing dependencies and even has smaller sizes:

julia> λ = 4//5
4//5

julia> n = 5
5

julia> Base.summarysize(JordanBlock(λ, n))
24

julia> using FillArrays

julia> Base.summarysize(Bidiagonal(Fill(λ,n), Fill(one(λ), n-1), :U))
56

[dlfivefifty]

But then you need to overload bandwidths etc.

[MikaelSlevinsky]

yes, i guess bandedmatrices is a dependency here?

[dlfivefifty]

I don't know what "dependencies" you are worried about. FillArrays.jl is already a dependency...

[dlfivefifty]

or at least is a dependency of ArrayLayouts.jl so effectively adding it as a dependency is a no-op

[MikaelSlevinsky]

ah ok

[MikaelSlevinsky]

Assuming we used const JordanBlock{T} = Bidiagonal{T, FillVector{T}}, aren't methods written here for JordanBlocks type piracy?

[MikaelSlevinsky]

Like exp(J::JordanBlock{T}) where T would modify Bidiagonal behaviour depending on whether or not this package is loaded.

[dlfivefifty]

It means you should overload exp in FillArrays.jl

none of the col/row support overloads are needed

[MikaelSlevinsky]

I guess exp(B::Bidiagonal{T, FillVector{T}} could live in FillArrays, but even then if you wanted to spell out the Toeplitz structure, you'd need to load ToeplitzMatrices in FillArrays. Maybe this is handled with an extension?

[dlfivefifty]

At the moment the only solution to this is to do type-piracy in ToeplitzMayrices.jl. But this type piracy (where a downstream package changes the return type to add structure) happens all over the place so is acceptable.