### A Pluto.jl notebook ###
# v0.14.4

using Markdown
using InteractiveUtils

# ╔═╡ 12bae146-12e7-4a78-afac-f5c2d6b86b66
begin
	using PlutoUI
	PlutoUI.TableOfContents(aside=true)
end

# ╔═╡ 7be5652e-3e5a-4199-a099-40f2da28053c
using LinearAlgebra, Arpack, LinearMaps, SparseArrays

# ╔═╡ 4b63c1d9-f043-448f-9e05-911f52d4227d
begin
	using Random
	Random.seed!(421)
	n=100
	A=Matrix(Symmetric(rand(n,n)))
	# Or: A = rand(5,5) |> t -> t + t'
	x=rand(n)
	k=10
end

# ╔═╡ f647d0dd-1fe4-42cf-b55c-38baa12f2db8
md"""
# Symmetric Eigenvalue Decomposition - Lanczos Method


If the matrix $A$ is large and sparse and/or if only some eigenvalues and their eigenvectors are desired, iterative methods are the methods of choice. For example, the power method can be useful to compute the eigenvalue with the largest modulus. The basic operation in the power method is matrix-vector multiplication, and this can be performed very fast if $A$ is sparse. Moreover, $A$ need not be stored in the computer -- the input for the algorithm can be just a function which, given some vector $x$, returns the product $Ax$.

An _improved_ version of the power method, which efficiently computes some eigenvalues (either largest in modulus or near some target value $\mu$) and the corresponding eigenvectors, is the Lanczos method.

For more details, see [I. Slapničar, Symmetric Matrix Eigenvalue Techniques, pp. 55.1-55.25](https://www.routledge.com/Handbook-of-Linear-Algebra/Hogben/p/book/9781138199897) and the references therein.

__Prerequisites__

The reader should be familiar with concepts of eigenvalues and eigenvectors, related perturbation theory, and algorithms. 

 
__Competences__

The reader should be able to recognise matrices which warrant use uf Lanczos method, to apply the method and to assess the accuracy of the solution.
"""

# ╔═╡ d5f270bd-94c1-4da6-a8e6-a53337488020
md"""
# Lanczos method

 $A$ is a real symmetric matrix of order $n$.

## Definitions

Given a nonzero vector $x$ and an index $k<n$, the __Krylov matrix__ is defined as

$$
K_k=\begin{bmatrix} x & Ax & A^2 x &\cdots & A^{k-1}x \end{bmatrix}.$$

__Krilov subspace__ is the subspace spanned by the columns of $K_k$.

## Facts

1. The Lanczos method is based on the following observation. If $K_k=XR$ is the $QR$ factorization of the matrix $K_k$, then the $k\times k$ matrix $T=X^T A X$ is tridiagonal. The matrices $X$ and $T$ can be computed by using only matrix-vector products in $O(kn)$ operations.

2. Let $T=Q\Lambda Q^T$ be the EVD of $T$. Then $\lambda_i$ approximate well some of the largest and smallest eigenvalues of $A$, and the columns of the matrix $U=XQ$ approximate the corresponding eigenvectors.

3. As $k$ increases, the largest (smallest) eigenvalues of the matrix $T_{1:k,1:k}$ converge towards some of the largest (smallest) eigenvalues of $A$ (due to the Cauchy interlace property). The algorithm can be redesigned to compute only largest or smallest eigenvalues. Also, by using shift and invert strategy, the method can be used to compute eigenvalues near some specified value. In order to obtain better approximations, $k$ should be greater than the number of required eigenvalues. On the other side, in order to obtain better accuracy and efficacy, $k$ should be as small as possible.

4. The last computed element, $\mu=T_{k+1,k}$, provides information about accuracy:

$$
\begin{aligned}
\|AU-U\Lambda\|_2&=\mu, \\
\|AU_{:,i}-\lambda_i U_{:,i}\|_2&=\mu |Q_{ki}|, \quad  i=1,\ldots,k.
\end{aligned}$$

Further, there are $k$ eigenvalues $\tilde\lambda_1,\ldots,\tilde\lambda_k$ of $A$ such that $|\lambda_i-\tilde\lambda_i|\leq \mu$, and for the corresponding eigenvectors, we have 

$$\sin2\Theta(U_{:,i},\tilde U_{:,i}) \leq \frac{2\mu}{\min_{j\neq i} |\lambda_i-\tilde \lambda_j|}.$$ 

5. In practical implementations, $\mu$ is usually used to determine the index $k$. 

6. The Lanczos method has inherent numerical instability in the floating-point arithmetic: since the Krylov vectors are, in fact, generated by the power method, they converge towards an eigenvector of $A$. Thus, as $k$ increases, the Krylov vectors become more and more parallel, and the recursion in the function `Lanczos()` becomes numerically unstable and the computed columns of $X$ cease to be sufficiently orthogonal. This affects both the convergence and the accuracy of the algorithm. For example, several eigenvalues of $T$ may converge towards a simple eigenvalue of $A$ (the, so called, __ghost eigenvalues__).

7. The loss of orthogonality is dealt with by using the __full reorthogonalization__ procedure: in each step, the new ${\bf z}$ is orthogonalized against all previous columns of $X$, that is, in function `Lanczos()`, the formula 
```
z=z-Tr.dv[i]*X[:,i]-Tr.ev[i-1]*X[:,i-1]
```
is replaced by
```
z=z-sum(dot(z,Tr.dv[i])*X[:,i]-Tr.ev[i-1]*X[:,i-1]
```
  
To obtain better orthogonality, the latter formula is usually executed twice. The full reorthogonalization raises the operation count to $O(k^2n)$.

8. The __selective reorthogonalization__ is the procedure in which the current $z$ is orthogonalized against some selected columns of $X$, in order to attain sufficient numerical stability and not increase the operation count too much. The details are very subtle and can be found in the references.
  
9. The Lanczos method is usually used for sparse matrices. Sparse matrix $A$ is stored in the sparse format in which only values and indices of nonzero elements are stored. The number of operations required to multiply some vector by $A$ is also proportional to the number of nonzero elements.
  
10. The function `eigs()` implements Lanczos method real for symmetric matrices and more general Arnoldi method for general matrices.
"""

# ╔═╡ 97e0dbf8-b9be-4503-af5d-4cf6e86311eb
function Lanczos(A::Array{T}, x::Vector{T}, k::Int) where T
    n=size(A,1)
    X=Array{T}(undef,n,k)
    dv=Array{T}(undef,k)
    ev=Array{T}(undef,k-1)
    X[:,1]=x/norm(x)
    for i=1:k-1
        z=A*X[:,i]
        dv[i]=X[:,i]⋅z
        # Three-term recursion
        if i==1
            z=z-dv[i]*X[:,i]
        else
            # z=z-dv[i]*X[:,i]-ev[i-1]*X[:,i-1]
            # Full reorthogonalization - once or even twice
            z=z-sum([(z⋅X[:,j])*X[:,j] for j=1:i])
            z=z-sum([(z⋅X[:,j])*X[:,j] for j=1:i])
        end
        μ=norm(z)
        if μ==0
            Tr=SymTridiagonal(dv[1:i-1],ev[1:i-2])
            return eigvals(Tr), X[:,1:i-1]*eigvecs(Tr), X[:,1:i-1], μ
        else
            ev[i]=μ
            X[:,i+1]=z/μ
        end
    end
    # Last step
    z=A*X[:,end]
    dv[end]=X[:,end]⋅z
    z=z-dv[end]*X[:,end]-ev[end]*X[:,end-1]
    μ=norm(z)
    Tr=SymTridiagonal(dv,ev)
    eigvals(Tr), X*eigvecs(Tr), X, μ
end

# ╔═╡ 219ce78b-8bd3-4df3-93df-c08fad30e33f
λ,U,X,μ=Lanczos(A,x,10)

# ╔═╡ 1ef82905-e422-4644-ad00-26c448cb0e3a
# Orthogonality of X
norm(X'*X-I)

# ╔═╡ 595b22a5-4456-41fb-b4d6-861833bc6d47
# Tridiagonalization
X'*A*X

# ╔═╡ b2e762d4-650c-4846-ae30-32614b516fe3
# Residual
norm(A*U-U*Diagonal(λ)), μ

# ╔═╡ bd2b8d59-0525-4db5-99cd-e2f98f735eda
U'*A*U

# ╔═╡ 6791903f-0cf7-4bbc-adfb-ff9e80b373dc
# Orthogonality of U
norm(U'*U-I)

# ╔═╡ 218a2a21-e391-4a90-a96a-71bbb0d2f895
# Full eigenvalue decomposition
λeigen,Ueigen=eigen(A);

# ╔═╡ e23c1bc6-c4d8-4c5d-9a83-2b00c5d93b25
# ?eigs

# ╔═╡ fc00c272-d66f-42be-9964-7a934b32015c
# Lanczos method from Arpack.jl
λeigs,Ueigs=eigs(A; nev=k, which=:LM, ritzvec=true, v0=x)

# ╔═╡ 91eec1a8-413c-4960-bb52-8da2435e9e4a
[λ λeigs λeigen[sortperm(abs.(λeigen),rev=true)[1:k]] ]

# ╔═╡ aa20cf19-2bc8-4831-ae26-fb883614c63f
md"""
We see that `eigs()` computes `k` eigenvalues with largest modulus. What eigenvalues did `Lanczos()` compute?
"""

# ╔═╡ 9ce7a012-723d-4184-8ef3-00f468e61281
for i=1:k
    println(minimum(abs,λeigen.-λ[i]))
end

# ╔═╡ ed7d1400-e6d9-44e2-a633-7f6b3df74272
md"""
Conslusion is that the naive implementation of Lanczos is not enough. However, it is fine, when all eigenvalues are computed. Why?
"""

# ╔═╡ d2f31110-5587-4dd3-a5cf-cc3e046e1ee0
λall,Uall,Xall,μall=Lanczos(A,x,100)

# ╔═╡ 04763c96-3a2b-4092-8969-99c18cc8fd54
# Residual and relative errors 
norm(A*Uall-Uall*Diagonal(λall)), norm((λeigen-λall)./λeigen)

# ╔═╡ 95ef94d2-129f-4dbb-b705-9a2c8660e22e
md"""
# Operator version

We can use Lanczos method with operator which, given vector `x`, returns the product `A*x`. We use the function `LinearMap()` from the package [LinearMaps.jl](https://github.com/Jutho/LinearMaps.jl)
"""

# ╔═╡ b8716ffe-8405-4018-bdb9-29c8cd2243dc
# ?LinearMap

# ╔═╡ 3bd1283d-91fa-4fdd-bcf8-946ac83b848c
# Operator from the matrix
C=LinearMap(A)

# ╔═╡ 0f8596c1-ee32-4b0d-b13b-d36fb611c99e
begin
	λC,UC=eigs(C; nev=k, which=:LM, ritzvec=true, v0=x)
	λeigs-λC
end

# ╔═╡ 05cd7b4e-18cf-485c-8113-4855f22ac8a0
md"""
Here is an example of `LinearMap()` with the function. 
"""

# ╔═╡ a8aeb3ce-5a67-48c7-9c15-8cc11e7ec39b
f(x)=A*x

# ╔═╡ 86bcfb51-6508-4644-ab3d-ee5d9675e1d6
D=LinearMap(f,n,issymmetric=true)

# ╔═╡ 5c08501f-010a-4d17-96b7-e339b0299262
begin
	λD,UD=eigs(D, nev=k, which=:LM, ritzvec=true, v0=x)
	λeigs-λD
end

# ╔═╡ 8a8d59cb-0a69-412e-948e-8110f04419fe
md"""
# Sparse matrices
"""

# ╔═╡ f63b3be6-f5d4-49c7-b8c7-a03826f21ad4
# ?sprand

# ╔═╡ 7beaf848-ad66-47ff-9411-b4ea94476f38
# Generate a sparse symmetric matrix
C₁=sprand(n,n,0.05) |> t -> t+t'

# ╔═╡ d468bf54-54ad-4f07-87d7-ce3d666471aa
issymmetric(C₁)

# ╔═╡ 0f9a4d39-c00c-44b2-8831-eeca1b34e79c
eigs(C₁; nev=k, which=:LM, ritzvec=true, v0=x)

# ╔═╡ a2384a29-511e-43b5-831d-b5b9c2e85ef0


# ╔═╡ Cell order:
# ╟─12bae146-12e7-4a78-afac-f5c2d6b86b66
# ╟─f647d0dd-1fe4-42cf-b55c-38baa12f2db8
# ╟─d5f270bd-94c1-4da6-a8e6-a53337488020
# ╠═7be5652e-3e5a-4199-a099-40f2da28053c
# ╠═97e0dbf8-b9be-4503-af5d-4cf6e86311eb
# ╠═4b63c1d9-f043-448f-9e05-911f52d4227d
# ╠═219ce78b-8bd3-4df3-93df-c08fad30e33f
# ╠═1ef82905-e422-4644-ad00-26c448cb0e3a
# ╠═595b22a5-4456-41fb-b4d6-861833bc6d47
# ╠═b2e762d4-650c-4846-ae30-32614b516fe3
# ╠═bd2b8d59-0525-4db5-99cd-e2f98f735eda
# ╠═6791903f-0cf7-4bbc-adfb-ff9e80b373dc
# ╠═218a2a21-e391-4a90-a96a-71bbb0d2f895
# ╠═e23c1bc6-c4d8-4c5d-9a83-2b00c5d93b25
# ╠═fc00c272-d66f-42be-9964-7a934b32015c
# ╠═91eec1a8-413c-4960-bb52-8da2435e9e4a
# ╟─aa20cf19-2bc8-4831-ae26-fb883614c63f
# ╠═9ce7a012-723d-4184-8ef3-00f468e61281
# ╟─ed7d1400-e6d9-44e2-a633-7f6b3df74272
# ╠═d2f31110-5587-4dd3-a5cf-cc3e046e1ee0
# ╠═04763c96-3a2b-4092-8969-99c18cc8fd54
# ╟─95ef94d2-129f-4dbb-b705-9a2c8660e22e
# ╠═b8716ffe-8405-4018-bdb9-29c8cd2243dc
# ╠═3bd1283d-91fa-4fdd-bcf8-946ac83b848c
# ╠═0f8596c1-ee32-4b0d-b13b-d36fb611c99e
# ╟─05cd7b4e-18cf-485c-8113-4855f22ac8a0
# ╠═a8aeb3ce-5a67-48c7-9c15-8cc11e7ec39b
# ╠═86bcfb51-6508-4644-ab3d-ee5d9675e1d6
# ╠═5c08501f-010a-4d17-96b7-e339b0299262
# ╟─8a8d59cb-0a69-412e-948e-8110f04419fe
# ╠═f63b3be6-f5d4-49c7-b8c7-a03826f21ad4
# ╠═7beaf848-ad66-47ff-9411-b4ea94476f38
# ╠═d468bf54-54ad-4f07-87d7-ce3d666471aa
# ╠═0f9a4d39-c00c-44b2-8831-eeca1b34e79c
# ╠═a2384a29-511e-43b5-831d-b5b9c2e85ef0
