API

mutable struct KPMExpansion{T<:AbstractFloat, Tfft<:FFTW.r2rFFTWPlan}

A type to represent the Chebyshev polynomial expansion used in the KPM algorithm.

Fields

• M::Int: Order of the Chebyshev polynomial expansion.
• bounds::NTuple{2,T}: Bounds on eigenspectrum in the KPM algorithm.
• buf::Vector{T}: The first M elements of this vector of the Chebyshev expansion coefficients.
• r2rplan::Tfft: Plan for performing DCT to efficiently calculate the Chebyshev expansion coefficients via Chebyshev-Gauss quadrature.

KPMExpansion(func::Function, bounds, M::Int, N::Int = 2*M)

Initialize an instance of the KPMExpansion type to approximate the univariate function func, called as func(x), with a order M Chebyshev polynomial expansion on the interval bounds[1] < x bounds[2]. Here, N ≥ M is the number of points at which func is evaluated on that specified interval, which are then used to calculate the expansion coeffiencents via Chebyshev-Gauss quadrature.

kpm_update!(
    kpm_expansion::KPMExpansion{T}, func::Function, bounds, M::Int, N::Int = 2*M
) where {T<:AbstractFloat}

In-place update an instance of KPMExpansion to reflect new values for eigenspectrum bounds, expansion order M, and number of points at which the expanded function is evaluated when computing the expansion coefficients. This includes recomputing the expansion coefficients.

kpm_update_bounds!(
    kpm_expansion::KPMExpansion{T}, func::Function, bounds
) where {T<:AbstractFloat}

In-place update an instance of KPMExpansion to reflect new values for eigenspectrum bounds, recomputing the expansion coefficients.

kpm_update_order!(
    kpm_expansion::KPMExpansion, func::Function, M::Int, N::Int = 2*M
)

In-place update the expansion order M for an instance of KPMExpansion, recomputing the expansion coefficients. It is also possible to udpate the number of point N the function func is evaluated at to calculate the expansion coefficients.

kpm_coefs(func::Function, bounds, M::Int, N::Int = 2*M)

Calculate and return the Chebyshev expansion coefficients. Refer to kpm_coefs! for more information.

kpm_coefs!(
    coefs::AbstractVector{T}, func::Function, bounds,
    buf::AbstractVector{T} = zeros(T, 2*length(coefs)),
    r2rplan = FFTW.plan_r2r!(buf, FFTW.REDFT10)
) where {T<:AbstractFloat}

Calculate and record the Chebyshev polynomial expansion coefficients to order M in the vector ceofs for the function func on the interval (bounds[1], bounds[2]). Let length(buf) be the number of evenly spaced points on the interval for which func is evaluated when performing Chebyshev-Gauss quadrature to compute the Chebyshev polynomial expansion coefficients.

kpm_moments(
    M::Int, A, bounds, R::T,
    tmp = zeros(eltype(R), size(R)..., 3)
) where {T<:AbstractVecOrMat}

If $R$ is a vector, calculate and return the first $M$ moments as

\[\mu_m = \langle R | T_m(A^\prime) | R \rangle\]

where $A^\prime$ is the rescaled version of $A$ using the bounds. If $R$ is a matrix then calculate the moments as

\[\mu_m = \frac{1}{N} \sum_{n=1}^N \langle R_n | T_m(A^\prime) | R_n \rangle,\]

where $| R_n \rangle$ is the n'th column of $R$.

kpm_moments(
    M::Int, A, bounds, U::T, V::T,
    tmp = zeros(eltype(V), size(V)..., 3)
) where {T<:AbstractVecOrMat}

If $U$ and $V$ are vector then calculate and return the first $M$ moments

\[\mu_m = \langle U | T_m(A^\prime) | V \rangle,\]

where $A^\prime$ is the rescaled version of $A$ using the bounds. If $U$ and $V$ are matrices then calculate the moments as

\[\mu_m = \frac{1}{N} \sum_{n=1}^N \langle U_n | T_m(A^\prime) | V_n \rangle,\]

where $| U_n \rangle$ and $| V_n \rangle$ are are the n'th columns of each matrix.

kpm_moments(
    A, kpm_expansion::KPMExpansion,, R::T,
    tmp = zeros(eltype(R), size(R)..., 3)
) where {T<:AbstractVecOrMat}

If $R$ is a vector, calculate and return the first $kpm_expansion.M$ moments as

\[\mu_m = \langle R | T_m(A^\prime) | R \rangle\]

where $A^\prime$ is the rescaled version of $A$ using the kpm_expansion.bounds. If $R$ is a matrix then calculate the moments as

\[\mu_m = \frac{1}{N} \sum_{n=1}^N \langle R_n | T_m(A^\prime) | R_n \rangle,\]

where $| R_n \rangle$ is the n'th column of $R$.

kpm_moments(
    A, kpm_expansion::KPMExpansion, U::T, V::T,
    tmp = zeros(eltype(V), size(V)..., 3)
) where {T<:AbstractVecOrMat}

If $U$ and $V$ are vector then calculate and return the first kpm_expansion.M moments

\[\mu_m = \langle U | T_m(A^\prime) | V \rangle,\]

where $A^\prime$ is the rescaled version of $A$ using the kpm_expansion.bounds. If $U$ and $V$ are matrices then calculate the moments as

\[\mu_m = \frac{1}{N} \sum_{n=1}^N \langle U_n | T_m(A^\prime) | V_n \rangle,\]

where $| U_n \rangle$ and $| V_n \rangle$ are are the n'th columns of each matrix.

kpm_moments!(
    μ::AbstractVector, A, bounds, R::T,
    tmp = zeros(eltype(R), size(R)..., 3)
) where {T<:AbstractVecOrMat}

If $R$ is a vector, calculate and return the first $M$ moments as

\[\mu_m = \langle R | T_m(A^\prime) | R \rangle\]

where $A^\prime$ is the rescaled version of $A$ using the bounds. If $R$ is a matrix then calculate the moments as

\[\mu_m = \frac{1}{N} \sum_{n=1}^N \langle R_n | T_m(A^\prime) | R_n \rangle,\]

where $| R_n \rangle$ is the n'th column of $R$.

kpm_moments!(
    μ::AbstractVector, A, bounds, U::T, V::T,
    tmp = zeros(eltype(V), size(V)..., 3)
) where {T<:AbstractVecOrMat}

If $U$ and $V$ are vector then calculate and return the first $M$ moments

\[\mu_m = \langle U | T_m(A^\prime) | V \rangle,\]

where $A^\prime$ is the rescaled version of $A$ using the bounds. If $U$ and $V$ are matrices then calculate the moments as

\[\mu_m = \frac{1}{N} \sum_{n=1}^N \langle U_n | T_m(A^\prime) | V_n \rangle,\]

where $| U_n \rangle$ and $| V_n \rangle$ are are the n'th columns of each matrix.

kpm_moments!(
    μ::AbstractVector, A, kpm_expansion::KPMExpansion, R::T,
    tmp = zeros(eltype(R), size(R)..., 3)
) where {T<:AbstractVecOrMat}

If $R$ is a vector, calculate and return the first $M$ moments as

\[\mu_m = \langle R | T_m(A^\prime) | R \rangle\]

where $A^\prime$ is the rescaled version of $A$ using the kpm_expansion.bounds. If $R$ is a matrix then calculate the moments as

\[\mu_m = \frac{1}{N} \sum_{n=1}^N \langle R_n | T_m(A^\prime) | R_n \rangle,\]

where $| R_n \rangle$ is the n'th column of $R$.

kpm_moments!(
    μ::AbstractVector, A, kpm_expansion::KPMExpansion, U::T, V::T,
    tmp = zeros(eltype(V), size(V)..., 3)
) where {T<:AbstractVecOrMat}

If $U$ and $V$ are vector then calculate and return the first $M$ moments

\[\mu_m = \langle U | T_m(A^\prime) | V \rangle,\]

where $A^\prime$ is the rescaled version of $A$ using the kpm_expansion.bounds. If $U$ and $V$ are matrices then calculate the moments as

\[\mu_m = \frac{1}{N} \sum_{n=1}^N \langle U_n | T_m(A^\prime) | V_n \rangle,\]

where $| U_n \rangle$ and $| V_n \rangle$ are are the n'th columns of each matrix.

kpm_density(
    N::Int, μ::T, W
) where {T<:AbstractVector}

Given the KPM moments $\mu$, evaluate the corresponding spectral density $\rho(\epsilon)$ at $N$ energy points $\epsilon$, where $W$ is the spectral bandwidth or the bounds of the spectrum. This function returns both $\rho(\epsilon)$ and $\epsilon$.

kpm_density!(
    ρ::AbstractVector{T},
    ϵ::AbstractVector{T},
    μ::AbstractVector{T},
    bounds,
    r2rplan = FFTW.plan_r2r!(ρ, FFTW.REDFT01)
) where {T<:AbstractFloat}

Given the KPM moments $\mu$, calculate the spectral density $\rho(\epsilon)$ and corresponding energies $\epsilon$ where it is evaluated. Here bounds is bounds the eigenspecturm, such that the bandwidth is given by W = bounds[2] - bounds[1].

kpm_density!(
    ρ::AbstractVector{T},
    ϵ::AbstractVector{T},
    μ::AbstractVector{T},
    W::T,
    r2rplan = FFTW.plan_r2r!(ρ, FFTW.REDFT01),
) where {T<:AbstractFloat}

Given the KPM moments $\mu$, calculate the spectral density $\rho(\epsilon)$ and corresponding energies $\epsilon$ where it is evaluated. Here $W$ is the spectral bandwidth.

kpm_dot(
    coefs::T, moments::T
) where {T<:AbstractVector}

Calculate the inner product

\[\langle c | \mu \rangle = \sum_{m=1}^M c_m \cdot \mu_m,\]

where $c_m$ are the KPM expansion coefficients, and $\mu_m$ are the KPM moments.

kpm_dot(
    A, coefs::AbstractVector, bounds, R::T,
    tmp = zeros(eltype(R), size(R)..., 3)
) where {T<:AbstractVecOrMat}

If $R$ is a single vector, then calculate the inner product

\[\begin{align*} S & = \langle R | F(A^\prime) | R \rangle \\ & = \sum_{m=1}^M \langle R | c_m T_m(A^\prime) | R \rangle \end{align*},\]

where $A^\prime$ is the scaled version of $A$ using the bounds. If $R$ is a matrix, then calculate

\[\begin{align*} S & = \langle R | F(A^\prime) | R \rangle \\ & = \frac{1}{N} \sum_{n=1}^N \sum_{m=1}^M \langle R_n | c_m T_m(A^\prime) | R_n \rangle \end{align*},\]

where $| R_n \rangle$ is a column of $R$.

kpm_dot(
    A, coefs::AbstractVector, bounds, U::T, V::T,
    tmp = zeros(eltype(V), size(V)..., 3)
) where {T<:AbstractVecOrMat}

If $U$ and $V$ are single vectors, then calculate the inner product

\[\begin{align*} S & = \langle U | F(A^\prime) | V \rangle \\ & = \sum_{m=1}^M \langle U | c_m T_m(A^\prime) | V \rangle \end{align*},\]

where $A^\prime$ is the scaled version of $A$ using the bounds. If $U$ and $V$ are matrices, then calculate

\[\begin{align*} S & = \langle U | F(A^\prime) | V \rangle \\ & = \frac{1}{N} \sum_{n=1}^N \sum_{m=1}^M \langle U_n | c_m T_m(A^\prime) | V_n \rangle \end{align*},\]

where $| U_n \rangle$ and $| V_n \rangle$ are the columns of each matrix.

kpm_dot(
    kpm_expansion::KPMExpansion, moments::AbstractVector
)

Calculate the inner product

\[\langle c | \mu \rangle = \sum_{m=1}^M c_m \cdot \mu_m,\]

where $c_m$ are the KPM expansion coefficients, and $\mu_m$ are the moments. Refer to kpm_moments to see how to calculate the moments.

kpm_dot(
    A, kpm_expansion::KPMExpansion, R::T,
    tmp = zeros(eltype(R), size(R)..., 3)
) where {T<:AbstractVecOrMat}

If $R$ is a single vector, then calculate the inner product

\[\begin{align*} S & = \langle R | F(A) | R \rangle \\ & = \sum_{m=1}^M \langle R | c_m T_m(A^\prime) | R \rangle \end{align*},\]

where $A^\prime$ is the scaled version of $A$ using the bounds. If $R$ is a matrix, then calculate

\[\begin{align*} S & = \langle R | F(A) | R \rangle \\ & = \frac{1}{N} \sum_{n=1}^N \sum_{m=1}^M \langle R_n | c_m T_m(A^\prime) | R_n \rangle \end{align*},\]

where $| R_n \rangle$ is a column of $R$.

kpm_dot(
    A, kpm_expansion::KPMExpansion, U::T, V::T,
    tmp = zeros(eltype(V), size(V)..., 3)
) where {T<:AbstractVecOrMat}

If $U$ and $V$ are single vectors, then calculate the inner product

\[\begin{align*} S & = \langle U | F(A) | V \rangle \\ & = \sum_{m=1}^M \langle U | c_m T_m(A^\prime) | V \rangle \end{align*},\]

where $A^\prime$ is the scaled version of $A$ using the bounds. If $U$ and $V$ are matrices, then calculate

\[\begin{align*} S & = \langle U | F(A) | V \rangle \\ & = \frac{1}{N} \sum_{n=1}^N \sum_{m=1}^M \langle U_n | c_m T_m(A^\prime) | V_n \rangle \end{align*},\]

where $| U_n \rangle$ and $| V_n \rangle$ are the columns of each matrix.

kpm_mul(
    A, coefs::AbstractVector, bounds, v::T, tmp = zeros(eltype(v), size(v)..., 3)
) where {T<:AbstractVecOrMat}

Evaluate and return the vector $v^\prime = F(A) \cdot v$ where $F(A)$ is represented by the Chebyshev expansion. For more information refer to kpm_mul!.

kpm_mul(A, kpm_expansion::KPMExpansion, v::T) where {T<:AbstractVecOrMat}

Evaluate and return the vector $v^\prime = F(A) \cdot v$ where $F(A)$ is represented by the Chebyshev expansion. For more information refer to kpm_mul!.

kpm_mul!(
    v′::T, A, coefs::AbstractVector, bounds, v::T,
    tmp = zeros(eltype(v), size(v)..., 3)
) where {T<:AbstractVecOrMat}

Evaluates $v^\prime = F(A) \cdot v$, writing the result to $v^\prime$, where $F(A)$ is represented by the Chebyshev expansion. Here A is either a function that can be called as A(u,v) to evaluate $u = A\cdot v$, modifying u in-place, or is a type for which the operation mul!(u, A, v) is defined. The vector coefs contains Chebyshev expansion coefficients to approximate $F(A)$, where the eigenspectrum of $A$ is contained in the interval (bounds[1], bounds[2]) specified by the bounds argument. The vector v is vector getting multiplied by the Chebyshev expansion for $F(A)$. Lastly, tmp is an array used to avoid dynamic memory allocations.

kpm_mul!(
    v′::T, A, kpm_expansion::KPMExpansion, v::T,
    tmp = zeros(eltype(v), size(v)..., 3)
) where {T<:AbstractVecOrMat}

Evaluates $v^\prime = F(A) \cdot v$, writing the result to $v^\prime$, where $F(A)$ is represented by the Chebyshev expansion. Here A is either a function that can be called as A(u,v) to evaluate $u = A\cdot v$, modifying u in-place, or is a type for which the operation mul!(u, A, v) is defined. Lastly, the array tmp is used to avoid dynamic memory allocations.

kpm_lmul!(
    A, coefs::AbstractVector, v::T, bounds,
    tmp = zeros(eltype(v), size(v)..., 3)
) where {T<:AbstractVecOrMat}

Evaluates $v = F(A) \cdot v$, modifying $v$ in-place, where $F(A)$ is represented by the Chebyshev expansion. Here A is either a function that can be called as A(u,v) to evaluate $u = A\cdot v$, modifying u in-place, or is a type for which the operation mul!(u, A, v) is defined. The vector coefs contains Chebyshev expansion coefficients to approximate $F(A)$, where the eigenspectrum of $A$ is contained in the interval (bounds[1], bounds[2]) specified by the bounds argument. The vector v is vector getting multiplied by the Chebyshev expansion for $F(A)$. Lastly, tmp is an array used to avoid dynamic memory allocations.

kpm_lmul!(
    A, kpm_expansion::KPMExpansion, v::T,
    tmp = zeros(eltype(v), size(v)..., 3)
) where {T<:AbstractVecOrMat}

Evaluates $v = F(A) \cdot v$, modifying $v$ in-place, where $F(A)$ is represented by the Chebyshev expansion. Here A is either a function that can be called as A(u,v) to evaluate $u = A\cdot v$, modifying u in-place, or is a type for which the operation mul!(u, A, v) is defined. Lastly, the array tmp is used to avoid dynamic memory allocations.

kpm_eval(x::AbstractFloat, coefs, bounds)

Evaluate $F(x)$ where $x$ is real number in the interval bounds[1] < x < bound[2], and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs.

kpm_eval(A::AbstractMatrix, coefs, bounds)

Evaluate and return the matrix $F(A),$ where $A$ is an operator with strictly real eigenvalues that fall in the interval (bounds[1], bounds[2]) specified by the bounds argument, and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs.

kpm_eval(x::T, kpm_expansion::KPMExpansion{T}) where {T<:AbstractFloat}

Evaluate $F(x)$ where $x$ is real number in the interval bounds[1] < x < bound[2], and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs.

kpm_eval(A::AbstractMatrix{T}, kpm_expansion::KPMExpansion{T}) where {T<:AbstractFloat}

Evaluate and return the matrix $F(A),$ where $A$ is an operator with strictly real eigenvalues and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs.

kpm_eval!(
    F::AbstractMatrix, A, coefs::AbstractVector, bounds,
    tmp = zeros(eltype(F), size(F)..., 4)
)

Evaluate and write the matrix $F(A)$ to F, where $A$ is an operator with strictly real eigenvalues that fall in the interval (bounds[1], bounds[2]) specified by the bounds argument, and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs. Lastly, tmp is used to avoid dynamic memory allocations.

kpm_eval!(
    F::AbstractMatrix, A, kpm_expansion::KPMExpansion,
    tmp = zeros(eltype(F), size(F)..., 3)
)

Evaluate and write the matrix $F(A)$ to F, where $A$ is an operator with strictly real eigenvalues and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs. Lastly, the array tmp is used to avoid dynamic memory allocations.

apply_jackson_kernel(coefs)

Return the Chebyshev expansion coefficients transformed by the Jackson kernel.

apply_jackson_kernel!(coefs)

Modify the Chebyshev expansion coefficients by applying the Jackson kernel to them.

apply_jackson_kernel!(kpm_expansion::KPMExpansion)

Modify the Chebyshev expansion coefficients by applying the Jackson kernel to them.

lanczos(niters, v, A, S = I, rng = Random.default_rng())

Use niters Lanczos iterations to find a truncated tridiagonal representation of $A\cdot S$, up to similarity transformation. Here, $A$ is any Hermitian matrix, while $S$ is both Hermitian and positive definite. Traditional Lanczos uses the identity matrix, $S = I$. The extension to non-identity matrices $S$ is as follows: Each matrix-vector product $A\cdot v$ becomes $(A S) \cdot v$, and each vector inner product $w^\dagger \cdot v$ becomes $w^\dagger \cdot S \cdot v$. The implementation below follows Wikipedia, and is the most stable of the four variants considered by Paige [1]. This implementation introduces additional vector storage so that each Lanczos iteration requires only one matrix-vector multiplication for $A$ and $S$, respectively.

This function returns a SymTridiagonal matrix. Note that the eigmin and eigmax routines have specialized implementations for a SymTridiagonal matrix type.

Similar generalizations of Lanczos have been considered in [2] and [3].

1. C. C. Paige, IMA J. Appl. Math., 373-381 (1972),

https://doi.org/10.1093%2Fimamat%2F10.3.373.

2. H. A. van der Vorst, Math. Comp. 39, 559-561 (1982),

https://doi.org/10.1090/s0025-5718-1982-0669648-0

3. M. Grüning, A. Marini, X. Gonze, Comput. Mater. Sci. 50, 2148-2156 (2011),

https://doi.org/10.1016/j.commatsci.2011.02.021.

lanczos!(
    αs::AbstractVector, βs::AbstractVector, v::AbstractVector,
    A, S = I;
    tmp::AbstractMatrix = zeros(eltype(v), length(v), 5),
    rng = Random.default_rng()
)

lanczos!(
    αs::AbstractMatrix, βs::AbstractMatrix, v::AbstractMatrix,
    A, S = I;
    tmp::AbstractArray = zeros(eltype(v), size(v)..., 5),
    rng = Random.default_rng()
)

Use Lanczos iterations to find a truncated tridiagonal representation of $A\cdot S$, up to similarity transformation. Here, $A$ is any Hermitian matrix, while $S$ is both Hermitian and positive definite. Traditional Lanczos uses the identity matrix, $S = I$. The extension to non-identity matrices $S$ is as follows: Each matrix-vector product $A\cdot v$ becomes $(A S) \cdot v$, and each vector inner product $w^\dagger \cdot v$ becomes $w^\dagger \cdot S \cdot v$. The implementation below follows Wikipedia, and is the most stable of the four variants considered by Paige [1]. This implementation introduces additional vector storage so that each Lanczos iteration requires only one matrix-vector multiplication for $A$ and $S$, respectively.

The number of Lanczos iterations performed equals niters = length(αs), and niters - 1 == length(βs). This function returns a SymTridiagonal matrix based on the contents of the vectors αs and βs. Note that the eigmin and eigmax routines have specialized implementations for a SymTridiagonal matrix type.

Note that if αs, βs and v are all matrices, then each column of v is treated as a seperate vector, and a vector of SymTridiagonal of length size(v,2) will be returned.

Similar generalizations of Lanczos have been considered in [2] and [3].

1. C. C. Paige, IMA J. Appl. Math., 373-381 (1972), https://doi.org/10.1093%2Fimamat%2F10.3.373.

2. H. A. van der Vorst, Math. Comp. 39, 559-561 (1982), https://doi.org/10.1090/s0025-5718-1982-0669648-0

3. M. Grüning, A. Marini, X. Gonze, Comput. Mater. Sci. 50, 2148-2156 (2011), https://doi.org/10.1016/j.commatsci.2011.02.021.