Notebook

The Original EigenDecompression.EigenDecompose¶

See

In [1]:

using LinearAlgebra
M = rand(100,100)#対角化したい行列
E, P = eigen(M)

Out[1]:

Eigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}
values:
100-element Vector{ComplexF64}:
 -2.6713385149783693 - 0.5238483074100407im
 -2.6713385149783693 + 0.5238483074100407im
  -2.380694749857707 - 1.5819003151310398im
  -2.380694749857707 + 1.5819003151310398im
 -2.3445049148115737 - 0.8348768899082996im
 -2.3445049148115737 + 0.8348768899082996im
 -2.1226846073905064 - 0.22379425662578242im
 -2.1226846073905064 + 0.22379425662578242im
  -2.049294063421365 - 1.6320269189857641im
  -2.049294063421365 + 1.6320269189857641im
 -2.0178436003856652 + 0.0im
 -1.9229326770807857 - 1.7886924334849146im
 -1.9229326770807857 + 1.7886924334849146im
                     ⋮
  1.9543269536275725 + 0.5040016799489412im
  1.9954076975615094 - 1.0426916311345158im
  1.9954076975615094 + 1.0426916311345158im
  2.0560135421866663 - 0.48888998497523045im
  2.0560135421866663 + 0.48888998497523045im
   2.416364338786246 + 0.0im
   2.486887950907518 - 1.224341365689113im
   2.486887950907518 + 1.224341365689113im
  2.6754719890150778 - 0.39289719212777263im
  2.6754719890150778 + 0.39289719212777263im
  3.0130359977740344 + 0.0im
    50.0991251975706 + 0.0im
vectors:
100×100 Matrix{ComplexF64}:
    0.0928227+0.0727977im   …  0.00480172+0.0im  -0.0897557+0.0im
   -0.0220687+0.0949872im      -0.0683635+0.0im   -0.107816+0.0im
    -0.062205-0.0221392im        0.150972+0.0im   -0.104096+0.0im
 -0.000929597-0.0558679im         0.21009+0.0im   -0.105664+0.0im
   -0.0858376+0.0687551im       0.0157017+0.0im  -0.0975072+0.0im
    0.0620079-0.0456313im   …    0.120079+0.0im   -0.104935+0.0im
    0.0910325-0.0859799im       0.0523112+0.0im  -0.0967437+0.0im
   -0.0549273+0.0307295im       0.0888397+0.0im   -0.106674+0.0im
     0.142685+0.0105723im       0.0982881+0.0im  -0.0831259+0.0im
    0.0213531-0.0756764im      -0.0218936+0.0im  -0.0993371+0.0im
   -0.0338946-0.108939im    …   0.0366303+0.0im  -0.0905221+0.0im
    0.0784134+0.105462im        0.0501874+0.0im  -0.0885164+0.0im
   -0.0326154+0.101677im        0.0967074+0.0im   -0.103145+0.0im
             ⋮              ⋱                    
    0.0293922-0.0230606im      -0.0531536+0.0im   -0.102552+0.0im
    0.0929652+0.0773336im        0.073754+0.0im  -0.0988448+0.0im
    0.0585327-0.0135158im   …   0.0340458+0.0im   -0.100852+0.0im
    0.0205471-0.0148881im       -0.205876+0.0im   -0.107353+0.0im
     0.174649-0.0368603im      -0.0668649+0.0im   -0.108233+0.0im
   -0.0567686-0.0222343im       0.0213794+0.0im   -0.100761+0.0im
      0.10441-0.0925521im       0.0107639+0.0im  -0.0950483+0.0im
    0.0780944+0.102672im    …   -0.166945+0.0im   -0.102682+0.0im
   -0.0634175-0.00904586im      0.0309103+0.0im   -0.102103+0.0im
     -0.11412-0.0948267im       0.0745764+0.0im  -0.0954268+0.0im
    0.0400251+0.0790165im       0.0617443+0.0im   -0.095069+0.0im
     0.068089-0.00059643im     -0.0551766+0.0im  -0.0948394+0.0im

In [2]:

exp(eigen(M))

MethodError: no method matching exp(::Eigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}})

Closest candidates are:
  exp(::Union{Float16, Float32, Float64})
   @ Base special\exp.jl:325
  exp(::Adjoint{T, <:AbstractMatrix} where T)
   @ LinearAlgebra D:\Julia-1.9.0-beta2\share\julia\stdlib\v1.9\LinearAlgebra\src\dense.jl:595
  exp(::Transpose{T, <:AbstractMatrix} where T)
   @ LinearAlgebra D:\Julia-1.9.0-beta2\share\julia\stdlib\v1.9\LinearAlgebra\src\dense.jl:596
  ...


Stacktrace:
 [1] top-level scope
   @ In[2]:1

In [3]:

module EigenDecompression

export EigenDecompose, eigDecomp
using LinearAlgebra
import Base.*, Base./

#対角化された行列型
struct EigenDecompose{T<:Number} <: AbstractMatrix{T}
    P::AbstractMatrix{T}
    D::Diagonal{T}
    invP::AbstractMatrix{T}
end

#普通のMatrixを対角化する
function eigDecomp(mat::AbstractMatrix)
    E, P = eigen(mat)
    EigenDecompose(P, Diagonal(E), inv(P))
end

#EigenDecompose型に対する関数
Base.exp(eig::EigenDecompose) = EigenDecompose(eig.P, exp(eig.D), eig.invP)
*(eig::EigenDecompose, vec::AbstractVector) = eig.P * eig.D * eig.invP * vec
*(eig::EigenDecompose, sc::Number) = EigenDecompose(eig.P, eig.D*sc, eig.invP)
/(eig::EigenDecompose, sc::Number) = EigenDecompose(eig.P, eig.D/sc, eig.invP)

#普通のMatrixに戻す
Base.Array(eig::EigenDecompose) = eig.P * eig.D * eig.invP

end

Out[3]:

Main.EigenDecompression

In [4]:

using .EigenDecompression

M = rand(100, 100)
eM = eigDecomp(M)
for i in 1:100
    v = rand(100)
    rnd = rand()
    @assert exp(M*rnd)*v ≈ exp(eM*rnd)*v
end

In [5]:

using BenchmarkTools

M = rand(100, 100);

#普通な方
function bench1(M)
    for i in 1:100
        v = rand(100)
        exp(M*rand())*v
    end
end

#今回実装した方
function bench2(M)
    eM = eigDecomp(M)
    for i in 1:100
        v = rand(100)
        exp(eM*rand())*v
    end
end

Out[5]:

bench2 (generic function with 1 method)

In [6]:

@benchmark bench1(M)

Out[6]:

BenchmarkTools.Trial: 58 samples with 1 evaluation.
 Range (min … max):  78.395 ms … 100.653 ms  ┊ GC (min … max): 2.43% … 3.65%
 Time  (median):     86.035 ms               ┊ GC (median):    2.82%
 Time  (mean ± σ):   87.298 ms ±   4.854 ms  ┊ GC (mean ± σ):  3.62% ± 1.60%

             ▂      █                                           
  ▄▄▁▁▄▄▁▁▁▁▁█▁▁█▄▆▁███▄▆▄▄▁██▁▄▄▆▄▁▄▁▁▆▁▄▄▆▁▁▄▆▁▁▁▁▁▁▄▆▁▁▁▁▁▄ ▁
  78.4 ms         Histogram: frequency by time         98.7 ms <

 Memory estimate: 53.80 MiB, allocs estimate: 1900.

In [7]:

@benchmark bench2(M)

Out[7]:

BenchmarkTools.Trial: 690 samples with 1 evaluation.
 Range (min … max):  6.103 ms … 15.411 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     6.824 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   7.234 ms ±  1.191 ms  ┊ GC (mean ± σ):  0.48% ± 2.20%

     █▅▂▁                                                     
  ▃▅▇████▇█▇▆▇▄▃▃▃▃▃▃▃▃▃▃▃▂▃▂▃▂▃▂▂▂▂▂▃▂▂▁▂▂▂▂▂▂▁▂▂▁▂▂▂▂▁▁▁▁▂ ▃
  6.1 ms         Histogram: frequency by time        12.3 ms <

 Memory estimate: 1.82 MiB, allocs estimate: 1728.

EigenDecomposedMatrices.EigenDecomposed¶

In [8]:

using LinearAlgebra
using BenchmarkTools

In [9]:

module EigenDecomposedMatrices

export EigenDecomposed

using LinearAlgebra
using Memoization

struct EigenDecomposed{
        T,
        TE<:AbstractVector{T},
        TP<:AbstractMatrix{T},
        TinvP<:AbstractMatrix{T}
    } <: AbstractMatrix{T}
    E::TE
    P::TP
    invP::TinvP
end

function EigenDecomposed(A::AbstractMatrix)
    E, P = eigen(A)
    invP = ishermitian(A) ? P' : inv(P)
    EigenDecomposed(E, P, invP)
end

LinearAlgebra.eigvals(ed::EigenDecomposed) = ed.E
LinearAlgebra.eigvecs(ed::EigenDecomposed) = ed.P
inveigvecs(ed::EigenDecomposed) = ed.invP
@memoize Base.parent(ed::EigenDecomposed) = eigvecs(ed) * Diagonal(eigvals(ed)) * inveigvecs(ed)
Base.convert(::Type{Array}, ed::EigenDecomposed) = convert(Array, parent(ed))
for op in (:eltype, :size)
    @eval Base.$op(ed::EigenDecomposed) = $op(eigvecs(ed))
end
Base.getindex(ed::EigenDecomposed, I...) = getindex(parent(ed), I...)

Base.:*(c::Number, ed::EigenDecomposed) = EigenDecomposed(c*eigvals(ed), eigvecs(ed), inveigvecs(ed))
Base.:*(ed::EigenDecomposed, c::Number) = EigenDecomposed(eigvals(ed)*c, eigvecs(ed), inveigvecs(ed))
Base.:\(c::Number, ed::EigenDecomposed) = EigenDecomposed(c\eigvals(ed), eigvecs(ed), inveigvecs(ed))
Base.:/(ed::EigenDecomposed, c::Number) = EigenDecomposed(eigvals(ed)/c, eigvecs(ed), inveigvecs(ed))
for T in (AbstractVector, AbstractMatrix)
    @eval function Base.:*(ed::EigenDecomposed, v::$T)
        E, P, invP = eigvals(ed), eigvecs(ed), inveigvecs(ed)
        P * (Diagonal(E) * (invP * v))
    end
end

function exp_old(ed::EigenDecomposed)
    E, P, invP = eigvals(ed), eigvecs(ed), inveigvecs(ed)
    expE = exp.(E)
    expA = P * Diagonal(expE) * invP 
    EigenDecomposed(expE, P, invP)
end

LinearAlgebra.lmul!(c::Number, ed::EigenDecomposed) = lmul!(c, eigvals(ed))
LinearAlgebra.rmul!(ed::EigenDecomposed, c::Number) = rmul!(eigvals(ed), c)
LinearAlgebra.ldiv!(c::Number, ed::EigenDecomposed) = ldiv!(c, eigvals(ed))
LinearAlgebra.rdiv!(ed::EigenDecomposed, c::Number) = rdiv!(eigvals(ed), c)

for op in (:exp, :log, :sin, :cos)
    opE = Symbol(op, "E")
    op_eigendecomposed = Symbol(op, "_eigendecomposed")
    op_eigendecomposed! = Symbol(op_eigendecomposed, "!")
    op_eigendecomposed!_doc =
        """
        $op_eigendecomposed!(Y, ed::EigenDecomposed, $opE=similar(ed.E), tmpY=similar(Y))

        returns the `$op` of `ed` and stores the result in `Y`, overwriting the existing value of `Y`. 
        It does not overwrite `ed` and uses `$opE` and `tmpY` as workspaces.
        """
    @eval begin
        @doc $op_eigendecomposed!_doc
        function $op_eigendecomposed!(Y, ed::EigenDecomposed, $opE=similar(ed.E), tmpY=similar(Y))
            E, P, invP = eigvals(ed), eigvecs(ed), inveigvecs(ed)
            @. $opE = $op.(E)
            mul!(tmpY, P, Diagonal($opE))
            mul!(Y, tmpY, invP)
        end
        $op_eigendecomposed(ed::EigenDecomposed) = $op_eigendecomposed!(similar(eigvecs(ed)), ed)
        Base.$op(ed::EigenDecomposed) = $op_eigendecomposed(ed)
    end
end

end

Out[9]:

Main.EigenDecomposedMatrices

In [10]:

?EigenDecomposedMatrices.exp_eigendecomposed!

Out[10]:

exp_eigendecomposed!(Y, ed::EigenDecomposed, expE=similar(ed.E), tmpY=similar(Y))

returns the exp of ed and stores the result in Y, overwriting the existing value of Y. It does not overwrite ed and uses expE and tmpY as workspaces.

In [11]:

?EigenDecomposedMatrices.log_eigendecomposed!

Out[11]:

log_eigendecomposed!(Y, ed::EigenDecomposed, logE=similar(ed.E), tmpY=similar(Y))

returns the log of ed and stores the result in Y, overwriting the existing value of Y. It does not overwrite ed and uses logE and tmpY as workspaces.

In [12]:

methods(EigenDecomposedMatrices.EigenDecomposed)

Out[12]:

# 2 methods for type constructor:

Main.EigenDecomposedMatrices.EigenDecomposed(E::TE, P::TP, invP::TinvP) where {T, TE<:AbstractVector{T}, TP<:AbstractMatrix{T}, TinvP<:AbstractMatrix{T}} in Main.EigenDecomposedMatrices at In[9]:14
Main.EigenDecomposedMatrices.EigenDecomposed(A::AbstractMatrix) in Main.EigenDecomposedMatrices at In[9]:19

In [13]:

methods(EigenDecomposedMatrices.EigenDecomposed{Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}})

Out[13]:

# 1 method for type constructor:

(var"#ctor-self#"::Type{Main.EigenDecomposedMatrices.EigenDecomposed{T, TE, TP, TinvP}} where {T, TE<:AbstractVector{T}, TP<:AbstractMatrix{T}, TinvP<:AbstractMatrix{T}})(E, P, invP) in Main.EigenDecomposedMatrices at In[9]:14

In [14]:

methodswith(EigenDecomposedMatrices.EigenDecomposed)

Out[14]:

21-element Vector{Method}:

*(ed::Main.EigenDecomposedMatrices.EigenDecomposed, c::Number) in Main.EigenDecomposedMatrices at In[9]:36
*(ed::Main.EigenDecomposedMatrices.EigenDecomposed, v::AbstractVector) in Main.EigenDecomposedMatrices at In[9]:40
*(ed::Main.EigenDecomposedMatrices.EigenDecomposed, v::AbstractMatrix) in Main.EigenDecomposedMatrices at In[9]:40
*(c::Number, ed::Main.EigenDecomposedMatrices.EigenDecomposed) in Main.EigenDecomposedMatrices at In[9]:35
/(ed::Main.EigenDecomposedMatrices.EigenDecomposed, c::Number) in Main.EigenDecomposedMatrices at In[9]:38
\(c::Number, ed::Main.EigenDecomposedMatrices.EigenDecomposed) in Main.EigenDecomposedMatrices at In[9]:37
convert(::Type{Array}, ed::Main.EigenDecomposedMatrices.EigenDecomposed) in Main.EigenDecomposedMatrices at In[9]:29
cos(ed::Main.EigenDecomposedMatrices.EigenDecomposed) in Main.EigenDecomposedMatrices at In[9]:78
eltype(ed::Main.EigenDecomposedMatrices.EigenDecomposed) in Main.EigenDecomposedMatrices at In[9]:31
exp(ed::Main.EigenDecomposedMatrices.EigenDecomposed) in Main.EigenDecomposedMatrices at In[9]:78
getindex(ed::Main.EigenDecomposedMatrices.EigenDecomposed, I...) in Main.EigenDecomposedMatrices at In[9]:33
log(ed::Main.EigenDecomposedMatrices.EigenDecomposed) in Main.EigenDecomposedMatrices at In[9]:78
parent(ed::Main.EigenDecomposedMatrices.EigenDecomposed) in Main.EigenDecomposedMatrices at D:\.julia\packages\Memoization\ut5GT\src\Memoization.jl:162
sin(ed::Main.EigenDecomposedMatrices.EigenDecomposed) in Main.EigenDecomposedMatrices at In[9]:78
size(ed::Main.EigenDecomposedMatrices.EigenDecomposed) in Main.EigenDecomposedMatrices at In[9]:31
eigvals(ed::Main.EigenDecomposedMatrices.EigenDecomposed) in Main.EigenDecomposedMatrices at In[9]:25
eigvecs(ed::Main.EigenDecomposedMatrices.EigenDecomposed) in Main.EigenDecomposedMatrices at In[9]:26
ldiv!(c::Number, ed::Main.EigenDecomposedMatrices.EigenDecomposed) in Main.EigenDecomposedMatrices at In[9]:55
lmul!(c::Number, ed::Main.EigenDecomposedMatrices.EigenDecomposed) in Main.EigenDecomposedMatrices at In[9]:53
rdiv!(ed::Main.EigenDecomposedMatrices.EigenDecomposed, c::Number) in Main.EigenDecomposedMatrices at In[9]:56
rmul!(ed::Main.EigenDecomposedMatrices.EigenDecomposed, c::Number) in Main.EigenDecomposedMatrices at In[9]:54

In [15]:

methods(EigenDecomposedMatrices.exp_eigendecomposed!)

Out[15]:

# 3 methods for generic function exp_eigendecomposed! from [36mMain.EigenDecomposedMatrices[39m:

exp_eigendecomposed!(Y, ed::Main.EigenDecomposedMatrices.EigenDecomposed) in Main.EigenDecomposedMatrices at In[9]:71
exp_eigendecomposed!(Y, ed::Main.EigenDecomposedMatrices.EigenDecomposed, expE) in Main.EigenDecomposedMatrices at In[9]:71
exp_eigendecomposed!(Y, ed::Main.EigenDecomposedMatrices.EigenDecomposed, expE, tmpY) in Main.EigenDecomposedMatrices at In[9]:71

In [16]:

methods(EigenDecomposedMatrices.log_eigendecomposed!)

Out[16]:

# 3 methods for generic function log_eigendecomposed! from [36mMain.EigenDecomposedMatrices[39m:

log_eigendecomposed!(Y, ed::Main.EigenDecomposedMatrices.EigenDecomposed) in Main.EigenDecomposedMatrices at In[9]:71
log_eigendecomposed!(Y, ed::Main.EigenDecomposedMatrices.EigenDecomposed, logE) in Main.EigenDecomposedMatrices at In[9]:71
log_eigendecomposed!(Y, ed::Main.EigenDecomposedMatrices.EigenDecomposed, logE, tmpY) in Main.EigenDecomposedMatrices at In[9]:71

In [17]:

A = [
    2 -1 0
    -1 2 -1
    0 -1 2
]

edA = EigenDecomposedMatrices.EigenDecomposed(A)

Out[17]:

3×3 Main.EigenDecomposedMatrices.EigenDecomposed{Float64, Vector{Float64}, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}:
  2.0          -1.0  -3.33067e-16
 -1.0           2.0  -1.0
 -3.33067e-16  -1.0   2.0

In [18]:

log(edA)

Out[18]:

3×3 Matrix{Float64}:
  0.51986   -0.623225  -0.173287
 -0.623225   0.346574  -0.623225
 -0.173287  -0.623225   0.51986

In [19]:

log(A)

Out[19]:

3×3 Matrix{Float64}:
  0.51986   -0.623225  -0.173287
 -0.623225   0.346574  -0.623225
 -0.173287  -0.623225   0.51986

In [20]:

log(edA) ≈ log(A)

Out[20]:

true

In [21]:

n = 2^8
M = 5I + randn(n, n)
v = randn(n)
c = randn()

edM = EigenDecomposedMatrices.EigenDecomposed(M)

Y = similar(eigvecs(edM))
expE = similar(eigvals(edM))
tmpY = similar(Y)

y = similar(eigvals(edM))
alpha = randn()
beta = randn();

In [22]:

edM

Out[22]:

256×256 Main.EigenDecomposedMatrices.EigenDecomposed{ComplexF64, Vector{ComplexF64}, Matrix{ComplexF64}, Matrix{ComplexF64}}:
    3.56334-1.65165e-13im  …  -0.0799542+7.35922e-14im
   -1.29471-3.33504e-14im        2.79704-1.73821e-15im
   0.188711-1.38904e-14im      -0.146065-7.04386e-16im
  -0.139691-1.2809e-13im        -0.11464+1.72234e-13im
   0.242261+2.90866e-14im       -1.44222+1.48799e-13im
   0.476316+3.58411e-14im  …    0.497164-2.26801e-13im
   0.550097-3.61273e-14im      -0.134351+1.41119e-13im
   0.717261+1.67825e-13im       0.895324-9.16857e-14im
   -1.04636-8.49634e-14im      -0.190769-2.40075e-13im
     1.6971-2.75795e-14im       0.770272-6.1418e-14im
    1.26714+5.92264e-14im  …   -0.567488-1.78777e-13im
 -0.0854576-1.88943e-14im      -0.822014+4.62595e-14im
   0.826429+6.51037e-14im       -1.53237-1.375e-13im
           ⋮               ⋱            ⋮
  -0.997608-2.67526e-14im      0.0412252+1.99661e-13im
 -0.0798052-1.49411e-14im  …    0.144672-1.3304e-13im
   -1.17006+1.4458e-14im       -0.431936+1.83672e-13im
 -0.0696381+8.73041e-14im      -0.825939-1.10264e-14im
   0.553973+1.47842e-14im        1.41854-6.53152e-14im
    1.05755+4.38409e-14im      -0.994764-9.63611e-14im
   0.143704+9.43705e-15im  …   -0.252374-1.1153e-13im
   0.089135-3.98563e-14im       -1.87436-1.5556e-14im
   -1.30653-7.54455e-14im        1.53995-1.05616e-13im
   0.642478-5.84129e-14im        1.36199-7.11765e-14im
    0.24546-2.13334e-14im       0.401485-2.4925e-14im
    -1.3641+7.4302e-14im   …     5.51035+9.54081e-14im

In [23]:

dump(edM)

Main.EigenDecomposedMatrices.EigenDecomposed{ComplexF64, Vector{ComplexF64}, Matrix{ComplexF64}, Matrix{ComplexF64}}
  E: Array{ComplexF64}((256,)) ComplexF64[-10.529011730788033 + 0.0im, -10.280415425445305 - 5.0635332402051345im, -10.280415425445305 + 5.0635332402051345im, -9.365531653942693 - 7.310921121421777im, -9.365531653942693 + 7.310921121421777im, -9.066423960276218 - 0.8998268797512045im, -9.066423960276218 + 0.8998268797512045im, -8.732971181835605 - 1.4609873383902336im, -8.732971181835605 + 1.4609873383902336im, -8.46124963643927 - 8.33273467746558im  …  18.632067916841073 - 4.578055223585328im, 18.632067916841073 + 4.578055223585328im, 18.93064239638064 - 7.073434597909069im, 18.93064239638064 + 7.073434597909069im, 19.20170503001391 - 3.658032128271094im, 19.20170503001391 + 3.658032128271094im, 20.930683959675996 - 0.7030285754956695im, 20.930683959675996 + 0.7030285754956695im, 21.184754820465372 - 3.0347453281790697im, 21.184754820465372 + 3.0347453281790697im]
  P: Array{ComplexF64}((256, 256)) ComplexF64[0.016739760214174452 + 0.0im 0.0479671595323605 - 0.01662350156735747im … 0.08231100940410303 + 0.03720618909538042im 0.08231100940410303 - 0.03720618909538042im; -0.05678105463615567 + 0.0im -0.02753855144836508 - 0.03255712650078475im … 0.01670703557346129 + 0.024130937819675107im 0.01670703557346129 - 0.024130937819675107im; … ; 0.023221059785412463 + 0.0im -0.017275722649348676 - 0.009015532372298877im … 0.049235288759355164 + 0.03697492625403338im 0.049235288759355164 - 0.03697492625403338im; -0.0578255784977335 + 0.0im -0.002679638667629752 - 0.007116991095201977im … 0.02694655853750013 + 0.0004816998404853853im 0.02694655853750013 - 0.0004816998404853853im]
  invP: Array{ComplexF64}((256, 256)) ComplexF64[0.3043258839874642 + 4.5105492045176234e-15im 0.11647828660858314 + 8.454429333244175e-15im … 0.15219311634828636 + 1.4654943925052066e-14im -0.17563710572420543 + 2.1649348980190553e-15im; 0.08912104015473599 + 0.06654733900857696im -0.08069004242424725 + 0.18453572669503732im … -0.2156137951599397 + 0.002174603423839006im 0.02083977906498702 - 0.0641592478208322im; … ; 0.0385888730719014 - 0.09966379283416266im 0.11941897834644964 + 0.033343849341449926im … -0.10247759199925907 + 0.07210367195681458im 0.048176002794752926 - 0.18025959146290277im; 0.03858887307190356 + 0.09966379283416041im 0.1194189783464544 - 0.03334384934145064im … -0.1024775919992489 - 0.07210367195681436im 0.048176002794748166 + 0.18025959146290615im]

In [24]:

M ≈ parent(edM) == Matrix(edM)

Out[24]:

true

In [25]:

M ≈ edM

Out[25]:

true

In [26]:

c*M ≈ c*edM ≈ edM*c

Out[26]:

true

In [27]:

c\M ≈ c\edM ≈ edM/c

Out[27]:

true

In [28]:

(
    exp(M)
    ≈ exp(edM)
    ≈ EigenDecomposedMatrices.exp_old(edM)
    ≈ EigenDecomposedMatrices.exp_eigendecomposed!(Y, edM)
    ≈ EigenDecomposedMatrices.exp_eigendecomposed!(Y, edM, expE, tmpY)
)

Out[28]:

true

In [29]:

@show typeof(y)
(
    exp(M) * v
    ≈ exp(edM) * v
    ≈ EigenDecomposedMatrices.exp_old(edM) * v
    ≈ mul!(y, EigenDecomposedMatrices.exp_eigendecomposed!(Y, edM), v)
    ≈ mul!(y, EigenDecomposedMatrices.exp_eigendecomposed!(Y, edM, expE, tmpY), v)
)

typeof(y) = Vector{ComplexF64}

Out[29]:

true

In [30]:

@btime edM = EigenDecomposedMatrices.EigenDecomposed(M);

  40.167 ms (26 allocations: 3.53 MiB)

In [31]:

@btime exp($M) * $v
@btime exp($edM) * $v
@btime EigenDecomposedMatrices.exp_old($edM) * $v
@btime mul!($y, EigenDecomposedMatrices.exp_eigendecomposed!($Y, $edM), $v)
@btime mul!($y, EigenDecomposedMatrices.exp_eigendecomposed!($Y, $edM, $expE, $tmpY), $v);

  10.090 ms (16 allocations: 3.01 MiB)
  2.712 ms (7 allocations: 2.01 MiB)
  2.696 ms (9 allocations: 2.02 MiB)
  2.540 ms (3 allocations: 1.00 MiB)
  2.447 ms (0 allocations: 0 bytes)

In [32]:

n2 = 2^8
M2 = Symmetric(5I + randn(n2, n2))
v2 = randn(n)
c2 = randn()

edM2 = EigenDecomposedMatrices.EigenDecomposed(M2)

Y2 = similar(eigvecs(edM2))
expE2 = similar(eigvals(edM2))
tmpY2 = similar(Y2)

y2 = similar(eigvals(edM2))
alpha2 = randn()
beta2 = randn();

In [33]:

edM2

Out[33]:

256×256 Main.EigenDecomposedMatrices.EigenDecomposed{Float64, Vector{Float64}, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}:
  6.03077    -0.43245     0.599883    …  -0.872836      0.980907
 -0.43245     7.07582    -0.965762       -2.46577       1.27889
  0.599883   -0.965762    6.28325        -0.897865      0.902268
  0.554185   -1.61176    -1.88607         1.52239      -0.800269
 -0.0394612  -0.472386   -2.16329        -0.666488      1.60583
 -0.684983    0.673263    1.44873     …  -1.11132      -1.06482
  0.544451    1.8252      0.567267        0.568297      1.60087
 -0.60012    -1.6835     -0.332066        0.151455     -1.76098
 -1.32028    -1.71558    -0.503773       -0.276013      1.46591
 -0.619846    0.546905   -1.79987         0.501719      0.114942
 -1.96161     0.575018   -0.648195    …   1.16492      -1.17966
 -0.303545   -0.670053    0.967682       -0.000749713   0.373673
  1.33041    -0.0295437  -0.717992        0.155946      1.85579
  ⋮                                   ⋱                 ⋮
 -0.590561    0.333331    0.250405        0.90949      -1.53465
 -0.748456    1.29363     0.1356      …  -1.60434       0.683966
  1.53178     0.157997    0.00791121     -1.8972        0.444282
  0.269089   -0.167004    1.49271        -0.630714     -2.68816
  2.1885     -1.43587    -1.18202        -1.17287       1.17593
 -0.787258   -0.512251   -0.294453        2.0117        0.0495256
 -0.550615   -0.564627    1.8215      …   0.254222     -0.881607
  1.39278     0.444108    0.341001       -0.411732      0.491344
  0.960274   -0.737275    0.302913        1.16593      -0.539727
 -1.638      -0.849312    0.108656       -0.763322      1.67947
 -0.872836   -2.46577    -0.897865        6.45844       0.754344
  0.980907    1.27889     0.902268    …   0.754344      5.11205

In [34]:

dump(edM2)

Main.EigenDecomposedMatrices.EigenDecomposed{Float64, Vector{Float64}, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}
  E: Array{Float64}((256,)) [-26.281246714455474, -26.079893599626686, -25.32277485066867, -25.0801104625073, -23.954875700929986, -23.74285195174835, -23.168019607432065, -22.60925439400774, -22.454572564359246, -21.705343122758094  …  31.569737139516377, 32.20706026446964, 32.53684276189552, 33.124400141791504, 33.974663809226215, 34.44625078387939, 34.499597457606136, 34.84570025054353, 36.07654166938045, 36.899884804719484]
  P: Array{Float64}((256, 256)) [0.1449279306214216 0.005919536496859418 … 0.00196042283590344 0.005370199356810865; 0.0977684655406359 -0.04973701370055792 … -0.0878992312227967 -0.05778779773302181; … ; 0.04308743597446929 0.04324139883354045 … -0.019127423566981328 -0.04138098814577151; -0.13325192873325928 -0.028583547120122572 … 0.12277657998369695 -0.056353074728779456]
  invP: Adjoint{Float64, Matrix{Float64}}
    parent: Array{Float64}((256, 256)) [0.1449279306214216 0.005919536496859418 … 0.00196042283590344 0.005370199356810865; 0.0977684655406359 -0.04973701370055792 … -0.0878992312227967 -0.05778779773302181; … ; 0.04308743597446929 0.04324139883354045 … -0.019127423566981328 -0.04138098814577151; -0.13325192873325928 -0.028583547120122572 … 0.12277657998369695 -0.056353074728779456]

In [35]:

M2 ≈ parent(edM2) == Matrix(edM2)

Out[35]:

true

In [36]:

M2 ≈ edM2

Out[36]:

true

In [37]:

c2*M2 ≈ c2*edM2 ≈ edM2*c2

Out[37]:

true

In [38]:

c2\M2 ≈ c2\edM2 ≈ edM2/c2

Out[38]:

true

In [39]:

(
    exp(M2)
    ≈ exp(edM2)
    ≈ EigenDecomposedMatrices.exp_old(edM2)
    ≈ EigenDecomposedMatrices.exp_eigendecomposed!(Y2, edM2)
    ≈ EigenDecomposedMatrices.exp_eigendecomposed!(Y2, edM2, expE2, tmpY2)
)

Out[39]:

true

In [40]:

@show typeof(y2)
(
    exp(M2) * v2
    ≈ exp(edM2) * v2
    ≈ EigenDecomposedMatrices.exp_old(edM2) * v2
    ≈ mul!(y2, EigenDecomposedMatrices.exp_eigendecomposed!(Y2, edM2), v2)
    ≈ mul!(y2, EigenDecomposedMatrices.exp_eigendecomposed!(Y2, edM2, expE2, tmpY2), v2)
)

typeof(y2) = Vector{Float64}

Out[40]:

true

In [41]:

@btime edM2 = EigenDecomposedMatrices.EigenDecomposed(M2);

  6.472 ms (14 allocations: 1.59 MiB)

In [42]:

@btime exp($M2) * $v2
@btime exp($edM2) * $v2
@btime EigenDecomposedMatrices.exp_old($edM2) * $v2
@btime mul!($y2, EigenDecomposedMatrices.exp_eigendecomposed!($Y2, $edM2), $v2)
@btime mul!($y2, EigenDecomposedMatrices.exp_eigendecomposed!($Y2, $edM2, $expE2, $tmpY2), $v2);

  7.427 ms (19 allocations: 2.60 MiB)
  726.000 μs (6 allocations: 1.00 MiB)
  719.700 μs (8 allocations: 1.01 MiB)
  690.300 μs (3 allocations: 514.17 KiB)
  662.900 μs (0 allocations: 0 bytes)

In [ ]: