Usage

Quick Start

using UniqueKronecker

# Evaluate the Jacobian of the degree-2 unique monomial map at x = [3, 5]
x = [3.0, 5.0]
J = unique_kronecker_jacobian(x, 2)
# 3×2 Matrix{Float64}:
#  6.0  0.0
#  5.0  3.0
#  0.0  10.0

Offline Precomputation with `DiffMatCache`

For repeated Jacobian evaluations (e.g., in a time-stepping loop or OpInf pipeline), precompute the sparse differentiation matrices once and reuse them:

r = 5     # reduced dimension
pmax = 3  # maximum polynomial degree

# Precompute all N_d^{(i)} blocks for i = 1, ..., pmax
cache = DiffMatCache(r, pmax)

# Now evaluate Jacobians efficiently
x = randn(r)
J2 = unique_kronecker_jacobian(x, 2, cache)  # degree-2 Jacobian
J3 = unique_kronecker_jacobian(x, 3, cache)  # degree-3 Jacobian

In-Place Evaluation

To avoid allocations in tight loops, use the in-place variant:

using LinearAlgebra: binomial

r = 4; deg = 3
cache = DiffMatCache(r, deg)

qi = binomial(r + deg - 1, deg)
J = zeros(qi, r)   # preallocate output

for x in eachcol(snapshot_matrix)
    unique_kronecker_jacobian!(J, x, deg, cache)
    # ... use J ...
end

Batch Evaluation over Snapshot Matrices

Compute the Jacobian at every column of a snapshot matrix in one call:

r = 3; k = 100; deg = 2
X = randn(r, k)  # k snapshots of dimension r

# Returns a Vector{Matrix{Float64}} of length k
Js = unique_kronecker_jacobian(X, deg)

# With a precomputed cache
cache = DiffMatCache(r, deg)
Js = unique_kronecker_jacobian(X, deg, cache)

Building Differentiation Matrices Directly

Access the sparse differentiation matrix blocks $N_d^{(i)}$ for inspection or custom computations:

# Individual blocks: N_1^{(2)}, N_2^{(2)} for r=2, i=2
Nblocks = diffmat_blocks(2, 2)
Nblocks[1]  # N_1^{(2)} — sparse 3×2
Nblocks[2]  # N_2^{(2)} — sparse 3×2

# Full concatenated matrix: N^{(i)} = [N_1^{(i)}, ..., N_r^{(i)}]
N = diffmat(2, 2)  # sparse 3×4

Computing the Coupling Matrix

For the ManTA-OpInf augmented ROM, compute the coupling matrix $C(\hat{s}) = \sum_{i=2}^{p} G^{(i)} \nabla_{\hat{s}} p_i^u(\hat{s})$:

r = 4; p = 3; j = 6

# Precompute differentiation matrices
cache = DiffMatCache(r, p)

# G^{(i)} = Vj' * V^{(i)} — precomputed projection matrices
G = [randn(j, binomial(r + i - 1, i)) for i in 2:p]

# Coupling matrix at a single point
x = randn(r)
C = coupling_matrix(x, G, cache)  # j × r matrix

# Batch: coupling matrix at each column of a snapshot matrix
X = randn(r, 100)
Cs = coupling_matrix(X, G, cache)  # Vector of j × r matrices

# The mass matrix for the augmented ROM
M = I(r) + C' * C   # symmetric positive definite

Specialized Degree-2 and Degree-3 Jacobians

For maximum performance at low degrees, use the specialized routines that bypass sparse matrix–vector products entirely:

x = randn(5)
J2 = unique_kronecker_jacobian2(x)  # degree-2, direct loop
J3 = unique_kronecker_jacobian3(x)  # degree-3, direct loop

These produce identical results to unique_kronecker_jacobian(x, 2) and unique_kronecker_jacobian(x, 3), respectively, but avoid sparse mul! overhead.

Integration with the ManTA-OpInf Pipeline

Below is a sketch of how the Jacobian functions integrate into the full ManTA-OpInf Algorithm 1 (Formulation A):

using UniqueKronecker
using LinearAlgebra

# --- Phase 1–2: POD basis and polynomial manifold (assumed given) ---
# V   ∈ R^{n × r}        — POD basis
# Vbar ∈ R^{n × q}       — polynomial lifting operator [V^{(2)}, ..., V^{(p)}]
# Vj  ∈ R^{n × j}        — selected columns of Vbar
# Shat ∈ R^{r × k}       — reduced snapshots
# Sdothat ∈ R^{r × k}    — reduced time derivatives

# --- Phase 3: Offline precomputation for coupling matrix ---
r = size(V, 2)
cache = DiffMatCache(r, p)

# Precompute G^{(i)} = Vj' * V^{(i)}
G = Vector{Matrix{Float64}}(undef, p - 1)
col_offset = 0
for i in 2:p
    qi = binomial(r + i - 1, i)
    Vi = Vbar[:, col_offset+1 : col_offset+qi]
    G[i-1] = Vj' * Vi
    col_offset += qi
end

# Evaluate coupling matrices at each snapshot
Cs = coupling_matrix(Shat, G, cache)

# --- Phase 4: Operator inference (Formulation A) ---
# Assemble targets for the enriched block
Rbar = hcat([Cs[ℓ] * Sdothat[:, ℓ] for ℓ in 1:k]...)

# Solve the two decoupled OpInf problems using D, R, Rbar ...

API Reference

Differentiation Matrices

UniqueKronecker.diffmat_blocks — Function

diffmat_blocks(n::Int, i::Int) -> Vector{SparseMatrixCSC{Int,Int}}

Compute the n sparse differentiation matrices $N_d^{(i)} \in \mathbb{R}^{q_i \times q_{i-1}}$ for $d = 1, \dots, n$, as defined in Proposition 3.6.

Each matrix satisfies $[N_d^{(i)}]_{\alpha,\beta} = \alpha_d$ when $\beta = \alpha - e_d$ and $\alpha_d \ge 1$, and zero otherwise.

The column d of the Jacobian $\nabla_{\hat{s}} p_i^u(\hat{s})$ is then $N_d^{(i)} \, p_{i-1}^u(\hat{s})$.

Arguments

n::Int: dimension of the vector (number of variables).
i::Int: degree of the unique monomial map (must be ≥ 1).

Returns

Nblocks::Vector{SparseMatrixCSC{Int,Int}}: vector of length n.

Example

julia> Nblocks = diffmat_blocks(2, 2)
2-element Vector{SparseMatrixCSC{Int64, Int64}}:
 ...
julia> Nblocks[1]   # N_1^{(2)}
3×2 SparseMatrixCSC{Int64, Int64} with 2 stored entries:
 2  ⋅
 ⋅  1
 ⋅  ⋅
julia> Nblocks[2]   # N_2^{(2)}
3×2 SparseMatrixCSC{Int64, Int64} with 2 stored entries:
 ⋅  ⋅
 1  ⋅
 ⋅  2

source

UniqueKronecker.diffmat — Function

diffmat(n::Int, i::Int) -> N::SparseMatrixCSC{Int,Int}

Compute the full differentiation matrix $N^{(i)} = [N_1^{(i)}, \dots, N_n^{(i)}] \in \mathbb{R}^{q_i \times n \, q_{i-1}}$, the horizontal concatenation of all per-variable differentiation blocks.

The Jacobian of the unique monomial map satisfies (Proposition 3.6):

\[\nabla_{\hat{s}} p_i^u(\hat{s}) = N^{(i)} \bigl(I_n \otimes p_{i-1}^u(\hat{s})\bigr).\]

Arguments

n::Int: dimension of the vector.
i::Int: degree of the unique monomial map (must be ≥ 1).

Returns

N::SparseMatrixCSC{Int,Int}: sparse matrix of size $q_i \times n q_{i-1}$.

source

Cache

UniqueKronecker.DiffMatCache — Type

DiffMatCache(n::Int, pmax::Int)

Precompute and cache the sparse differentiation matrix blocks $N_d^{(i)}$ for $d = 1,\dots,n$ and $i = 1,\dots, p_{\max}$. These are constant matrices that depend only on the variable dimension n and the degree, so they can be assembled once offline and reused for every Jacobian evaluation.

Fields

n::Int: variable dimension.
pmax::Int: maximum polynomial degree cached.
blocks::Vector{Vector{SparseMatrixCSC{Int,Int}}}: blocks[i] is a length-n vector of the N_d^{(i)} blocks for degree i.

Example

julia> cache = DiffMatCache(3, 4)   # n=3 variables, up to degree 4
julia> cache.blocks[2]              # the three N_d^{(2)} blocks
julia> cache.blocks[2][1]           # N_1^{(2)}

source

Jacobian Evaluation

UniqueKronecker.unique_kronecker_jacobian — Function

unique_kronecker_jacobian(x::AbstractVector{T}, i::Int,
                          [cache::DiffMatCache]) -> J::Matrix{T}

Compute the Jacobian $\nabla_{\hat{s}} p_i^u(\hat{s})$ of the unique degree-i monomial map evaluated at x, returning a dense $q_i \times n$ matrix.

The computation uses the formula from Proposition 3.6:

\[\nabla_{\hat{s}} p_i^u(\hat{s}) = N^{(i)} \bigl(I_n \otimes p_{i-1}^u(\hat{s})\bigr),\]

evaluated column-by-column as $N_d^{(i)} \, p_{i-1}^u(\hat{s})$ for $d = 1,\dots,n$ to avoid forming the full Kronecker product.

Arguments

x::AbstractVector{T}: the point at which to evaluate (length n).
i::Int: degree of the monomial map (must be ≥ 1).
cache::DiffMatCache (optional): precomputed differentiation matrices. If omitted, they are computed on the fly.

Returns

J::Matrix{T}: the $q_i \times n$ Jacobian matrix.

Example

julia> x = [3.0, 5.0]
julia> J = unique_kronecker_jacobian(x, 2)
3×2 Matrix{Float64}:
 6.0  0.0
 5.0  3.0
 0.0  10.0

source

unique_kronecker_jacobian(X::AbstractMatrix{T}, i::Int,
                          [cache::DiffMatCache]) -> Vector{Matrix{T}}

Compute the Jacobian of the unique degree-i monomial map at each column of X, returning a vector of $q_i \times n$ Jacobian matrices.

Arguments

X::AbstractMatrix{T}: matrix whose columns are the evaluation points (size $n \times k$).
i::Int: degree of the monomial map.
cache::DiffMatCache (optional): precomputed differentiation matrices.

Returns

Js::Vector{Matrix{T}}: Js[ℓ] is the Jacobian at X[:, ℓ].

source

UniqueKronecker.unique_kronecker_jacobian! — Function

unique_kronecker_jacobian!(J::AbstractMatrix{T}, x::AbstractVector{T}, i::Int,
                           cache::DiffMatCache) -> J

In-place version of unique_kronecker_jacobian. Writes into the pre-allocated $q_i \times n$ matrix J.

Arguments

J::AbstractMatrix{T}: output matrix of size $q_i \times n$.
x::AbstractVector{T}: evaluation point (length n).
i::Int: degree of the monomial map.
cache::DiffMatCache: precomputed differentiation matrices.

Returns

J (modified in place).

source

UniqueKronecker.unique_kronecker_jacobian2 — Function

unique_kronecker_jacobian2(x::AbstractVector{T}) -> Matrix{T}

Specialized Jacobian of $p_2^u(\hat{s}) = [\hat{s}_1^2, \hat{s}_1\hat{s}_2, \dots]$ for degree 2, computed directly without differentiation matrices.

source

UniqueKronecker.unique_kronecker_jacobian3 — Function

unique_kronecker_jacobian3(x::AbstractVector{T}) -> Matrix{T}

Specialized Jacobian of $p_3^u(\hat{s})$ for degree 3, computed directly.

source

Coupling Matrix

UniqueKronecker.coupling_matrix — Function

coupling_matrix(x::AbstractVector{T}, G::Vector{<:AbstractMatrix},
                cache::DiffMatCache; pstart::Int=2) -> C::Matrix{T}

Compute the coupling matrix

\[C(\hat{s}) = \sum_{i=p_{\text{start}}}^{p} G^{(i)} \nabla_{\hat{s}} p_i^u(\hat{s}) \in \mathbb{R}^{j \times n}\]

where $G^{(i)} = V_j^\top V^{(i)} \in \mathbb{R}^{j \times q_i}$ are precomputed projection matrices and pstart (default 2) is the starting degree.

Arguments

x::AbstractVector{T}: reduced state $\hat{s}$ (length n).
G::Vector{<:AbstractMatrix}: G[i-pstart+1] is $G^{(i)}$ of size $j \times q_i$.
cache::DiffMatCache: precomputed differentiation matrices.
pstart::Int: first degree in the summation (default 2).

Returns

C::Matrix{T}: the $j \times n$ coupling matrix.

source

coupling_matrix(X::AbstractMatrix{T}, G::Vector{<:AbstractMatrix},
                cache::DiffMatCache; pstart::Int=2) -> Vector{Matrix{T}}

Batch version: compute the coupling matrix at each column of X.

source