HyperDualMatrixTools.jl Documentation

This package provides an overloaded factorize and \ that work with hyperdual-valued arrays.

It is essentially base on the hyper dual type defined by the HyperDualNumbers.jl package.

Motivation

The idea is that for a hyperdual-valued matrix $M = A + \varepsilon_1 B + \varepsilon_2 C + \varepsilon_1 \varepsilon_2 D$, its inverse is given by $M^{-1} = (I - \varepsilon_1 A^{-1} B - \varepsilon_2 A^{-1} C - \varepsilon_1\varepsilon_2 A^{-1} (D - B A^{-1} C - C A^{-1} B)) A^{-1}$. Therefore, only the inverse of $A$ is required to evaluate the inverse of $M$. This package should be useful for evaluation of second derivatives of functions that use \ (e.g., with iterative solvers).

How it works

HyperDualMatrixTools.jl makes available a HyperDualFactors type which contains the factors of $A$ (i.e., the output of factorize, e.g., $L$ and $U$, or $Q$ and $R$) and the non-real parts of $M$ (i.e., $B$, $C$, and $D$). HyperDualMatrixTools.jl overloads factorize so that for a hyperdual-valued matrix M, factorize(M) creates an instance of HyperDualFactors. Finally, HyperDualMatrixTools.jl also overloads \ to efficiently solve hyperdual-valued linear systems of the type $M x = y$ by using the default \ with the factors of $A$ only.

Usage

Create your hyperdual-valued matrix M:

n = 4
A, B, C, D = rand(n, n), randn(n, n), rand(n, n), randn(n, n)
M = A + ε₁ * B + ε₂ * C + ε₁ε₂ * D
typeof(M)

# output

Array{HyperDualNumbers.Hyper{Float64},2}

(The ε₁, ε₂, and ε₁ε₂ constants are provided by HyperDualMatrixTools.jl for convenience.)

Factorize M:

Mf = factorize(M)
typeof(Mf)

# output

HyperDualFactors

Apply \ to solve systems of the type M * x = y

y = rand(n, 4) * [1.0, ε₁, ε₂, ε₁ε₂]
x = Mf \ y
M * x ≈ y

# output

true

Functions

LinearAlgebra.factorize — Function

factorize(M::SparseMatrixCSC{<:Hyper,<:Int})

Invokes factorize on just the real part of M and stores it along with the dual parts into a HyperDualFactors object.

source

factorize(M::Array{<:Hyper,2})

Invokes factorize on just the real part of M and stores it along with the dual parts into a HyperDualFactors object.

source

Base.:\ — Function

\(M::HyperDualFactors, y::AbstractVecOrMat{Float64})

Backsubstitution for HyperDualFactors. See HyperDualFactors for details.

source

\(M::HyperDualFactors, y::AbstractVecOrMat{Hyper256})

Backsubstitution for HyperDualFactors. See HyperDualFactors for details.

source

\(Af::Factorization{Float64}, y::AbstractVecOrMat{Hyper256})

Backsubstitution for HyperDual-valued RHS.

source

\(M::AbstractArray{<:Hyper,2}, y::AbstractVecOrMat)

Backslash (factorization and backsubstitution) for Dual-valued matrix M.

source

New types

HyperDualMatrixTools.HyperDualFactors — Type

HyperDualFactors

Container type to work efficiently with backslash on hyperdual-valued sparse matrices.

factorize(M) will create an instance containing

Af = factorize(realpart.(M)) — the factors of the real part
B = ε₁part.(M) — the $\varepsilon_1$ part
C = ε₂part.(M) — the $\varepsilon_2$ part
D = ε₁ε₂part.(M) — the $\varepsilon_1\varepsilon_2$ part

for a hyperdual-valued matrix M.

This is because only the factors of the real part are needed when solving a linear system of the type $M x = b$ for a hyperdual-valued matrix $M = A + \varepsilon_1 B + \varepsilon_2 C + \varepsilon_1 \varepsilon_2 D$. In fact, the inverse of $M$ is given by $M^{-1} = (I - \varepsilon_1 A^{-1} B - \varepsilon_2 A^{-1} C - \varepsilon_1\varepsilon_2 A^{-1} (D - B A^{-1} C - C A^{-1} B)) A^{-1}$.

source