API

Base.iterate(it::UniqueCombinationIterator, comb::Vector{Int})

Iterate over all combinations with repetition of n elements from {1, 2, ..., n}.

⊘(x::AbstractVector{T}, y::AbstractVector{T}) where T

Unique Kronecker product operation

Arguments

• x::AbstractVector{T}: vector to perform the unique Kronecker product
• y::AbstractVector{T}: vector to perform the unique Kronecker product

Returns

• unique Kronecker product

⊘(x::AbstractArray{T}...) where {T<:Number}

Generalized Kronecker product operator for multiple vectors.

Arguments

• x::AbstractArray{T}...: one or more vectors to perform the unique Kronecker product

Returns

• unique Kronecker product of all vectors

⊛(args::AbstractArray...)

Circulant Kronecker product operator for multiple arguments.

Arguments

• args::AbstractArray...: Vectors or matrices to perform the circulant Kronecker product.

Returns

• result: The circulant Kronecker product.

H2Q(H::AbstractArray) → Q

Convert the quadratic H operator into the Q operator

Arguments

• H::AbstractArray: Quadratic matrix of dimensions (n x n^2)

Returns

• the Q quadratic matrix of 3-dim tensor

Q2H(Q::AbstractArray) → H

Convert the quadratic Q operator into the H operator. The Q matrix is a 3-dim tensor with dimensions (n x n x n). Thus,

\[\mathbf{Q} = \begin{bmatrix} \mathbf{Q}_1 \\ \mathbf{Q}_2 \\ \vdots \\ \mathbf{Q}_n \end{bmatrix} \quad \text{where }~~ \mathbf{Q}_i \in \mathbb{R}^{n \times n}\]

Arguments

• Q::AbstractArray: Quadratic matrix in the 3-dim tensor form with dimensions (n x n x n)

Returns

• the H quadratic matrix

circulant_kron_snapshot_matrix(Xmat::AbstractArray{T}...) where {T<:Number}

Compute the circulant Kronecker product of a set of matrices, where each matrix is a snapshot matrix.

Arguments

• Xmat::AbstractArray{T}...: Snapshot matrices to compute the circulant Kronecker product.
• redundant::Bool=true: If true, return the circulant Kronecker product as-is. If false, return the circulant Kronecker product with redundant rows removed.

Returns

• result: The circulant Kronecker product of the snapshot matrices.

circulant_kronecker(args::AbstractArray...)

Circulant Kronecker product operation for multiple Kronecker products.

Arguments

• args::AbstractArray...: Vectors or matrices to perform the circulant Kronecker product.

Returns

• result: The circulant Kronecker product.

commat(m::Int, n::Int) → K

Create commutation matrix K of dimension m x n Magnus and Neudecker [1].

Arguments

• m::Int: row dimension of the commutation matrix
• n::Int: column dimension of the commutation matrix

Returns

• K: commutation matrix

Example

julia> commat(2,2)
4×4 SparseMatrixCSC{Float64, Int64} with 4 stored entries:
 1.0   ⋅    ⋅    ⋅
  ⋅    ⋅   1.0   ⋅
  ⋅   1.0   ⋅    ⋅ 
  ⋅    ⋅    ⋅   1.0

commat(m::Int) → K

Dispatch for the commutation matrix of dimensions (m, m)

Arguments

• m::Int: row and column dimension of the commutation matrix

Returns

• K: commutation matrix

duplicate(A::AbstractArray)

Duplicate the redundant polynomial coefficients in the matrix A with a unique set of coefficients and return the matrix with redundant coefficients.

Arguments

• A::AbstractArray: A matrix
• p::Int: The order of the polynomial, e.g., p = 2 for x ⊗ x.

Returns

• matrix with redundant coefficients

Example

julia> n = 2; P = rand(n,n); P *= P'; p = vec(P)
4-element Vector{Float64}:
 0.5085988756090203
 0.7704767970682769
 0.7704767970682769
 1.310279680309927

julia> Q = rand(n,n); Q *= Q'; q = vec(Q)
4-element Vector{Float64}:
 0.40940214810208353
 0.2295272821417254
 0.2295272821417254
 0.25503767587483905

julia> A2 = [p'; q']
2×4 Matrix{Float64}:
 0.257622  0.44721   0.0202203  0.247649
 1.55077   0.871029  0.958499   0.650717

julia> D2 = dupmat(2,2)
4×3 SparseMatrixCSC{Int64, Int64} with 4 stored entries:
 1  ⋅  ⋅
 ⋅  1  ⋅
 ⋅  1  ⋅
 ⋅  ⋅  1

julia> A2 * D2
2×3 Matrix{Float64}:
 0.508599  1.54095   1.31028
 0.409402  0.459055  0.255038

julia> duplicate(A2 * D2, 2)
2×4 Matrix{Float64}:
 0.508599  1.54095   0.0  1.31028
 0.409402  0.459055  0.0  0.255038

duplicate_symmetric(A::AbstractArray, p::Int)

Duplicate the redundant polynomial coefficients in the matrix A with a unique set of coefficients and return the matrix with redundant coefficients which are duplicated symmetrically. This guarantees that the operator is symmetric. The difference from duplicate is that we use the elimination matrix Lp and the symmetric commutation matrix Sp to multiply the A matrix.

Arguments

• A::AbstractArray: A matrix
• p::Int: The order of the polynomial, e.g., p = 2 for x ⊗ x.

Returns

• matrix with redundant coefficients duplicated symmetrically

Example

julia> n = 2; P = rand(n,n); P *= P'; p = vec(P)
4-element Vector{Float64}:
 0.5085988756090203
 0.7704767970682769
 0.7704767970682769
 1.310279680309927

julia> Q = rand(n,n); Q *= Q'; q = vec(Q)
4-element Vector{Float64}:
 0.40940214810208353
 0.2295272821417254
 0.2295272821417254
 0.25503767587483905

julia> A2 = [p'; q']
2×4 Matrix{Float64}:
 0.257622  0.44721   0.0202203  0.247649
 1.55077   0.871029  0.958499   0.650717

julia> D2 = dupmat(2,2)
4×3 SparseMatrixCSC{Int64, Int64} with 4 stored entries:
 1  ⋅  ⋅
 ⋅  1  ⋅
 ⋅  1  ⋅
 ⋅  ⋅  1

julia> A2 * D2
2×3 Matrix{Float64}:
 0.508599  1.54095   1.31028
 0.409402  0.459055  0.255038

julia> duplicate_symmetric(A2 * D2, 2)
2×4 Matrix{Float64}:
 0.508599  0.770477  0.770477  1.31028
 0.409402  0.229527  0.229527  0.255038

dupmat(n::Int, p::Int) -> Dp::SparseMatrixCSC{Int}

Create a duplication matrix of order p for a vector of length n Magnus and Neudecker [1].

Arguments

• n::Int: The length of the vector.
• p::Int: The order of the duplication matrix, e.g., p = 2 for x ⊗ x.

Output

• Dp::SparseMatrixCSC{Int}: The duplication matrix of order p.

Example

julia> dupmat(2,2)
4×3 SparseMatrixCSC{Int64, Int64} with 4 stored entries:
 1  ⋅  ⋅
 ⋅  1  ⋅
 ⋅  1  ⋅
 ⋅  ⋅  1

elimat(n::Int, p::Int) -> Lp::SparseMatrixCSC{Int}

Create an elimination matrix of order p for a vector of length n Magnus and Neudecker [1].

Arguments

• n::Int: The length of the vector.
• p::Int: The order of the elimination matrix, e.g., p = 2 for x ⊗ x.

Output

• Lp::SparseMatrixCSC{Int}: The elimination matrix of order p.

Example

julia> elimat(2,2)
3×4 SparseMatrixCSC{Int64, Int64} with 3 stored entries:
 1  ⋅  ⋅  ⋅
 ⋅  1  ⋅  ⋅
 ⋅  ⋅  ⋅  1

eliminate(A::AbstractArray, p::Int)

Eliminate the redundant polynomial coefficients in the matrix A and return the matrix with unique coefficients.

Arguments

• A::AbstractArray: A matrix
• p::Int: The order of the polynomial, e.g., p = 2 for x ⊗ x.

Returns

• matrix with unique coefficients

Example

julia> n = 2; P = rand(n,n); P *= P'; p = vec(P)
4-element Vector{Float64}:
 0.5085988756090203
 0.7704767970682769
 0.7704767970682769
 1.310279680309927

julia> Q = rand(n,n); Q *= Q'; q = vec(Q)
4-element Vector{Float64}:
 0.40940214810208353
 0.2295272821417254
 0.2295272821417254
 0.25503767587483905

julia> A2 = [p'; q']
2×4 Matrix{Float64}:
 0.257622  0.44721   0.0202203  0.247649
 1.55077   0.871029  0.958499   0.650717

julia> eliminate(A2, 2)
2×3 Matrix{Float64}:
 0.508599  1.54095   1.31028
 0.409402  0.459055  0.255038

extractF(F::Union{SparseMatrixCSC,VecOrMat}, r::Int) → F

Extracting the F matrix for POD basis of dimensions (N, r)

Arguments

• F: F matrix
• r: reduced order

Returns

• extracted F matrix

extractH(H::Union{SparseMatrixCSC,VecOrMat}, r::Int) → H

Extracting the H matrix for POD basis of dimensions (N, r)

Arguments

• H: H matrix
• r: reduced order

Returns

• extracted H matrix

insertF(Fi::Union{SparseMatrixCSC,VecOrMat}, N::Int) → F

Inserting values into the F matrix for higher dimensions

Arguments

• Fi: F matrix to insert
• N: the larger order

Returns

• inserted F matrix

insertH(Hi::Union{SparseMatrixCSC,VecOrMat}, N::Int) → H

Inserting values into the H matrix for higher dimensions

Arguments

• Hi: H matrix to insert
• N: the larger order

Returns

• inserted H matrix

insert2bilin(X::Union{SparseMatrixCSC,VecOrMat}, N::Int, p::Int) → BL

Inserting the values into the bilinear matrix (N) for higher dimensions

Arguments

• X: bilinear matrix to insert
• N: the larger order

Returns

• Inserted bilinear matrix

insert2randF(Fi::Union{SparseMatrixCSC,VecOrMat}, N::Int) → F

Inserting values into the F matrix for higher dimensions

Arguments

• Fi: F matrix to insert
• N: the larger order

Returns

• inserted F matrix

invec(r::AbstractArray, m::Int, n::Int) → r

Inverse vectorization.

Arguments

• r::AbstractArray: the input vector
• m::Int: the row dimension
• n::Int: the column dimension

Returns

• the inverse vectorized matrix

kron_snapshot_matrix(Xmat::AbstractArray{T}, p::Int) where {T<:Number}

Take the p-order Kronecker product of each state of the snapshot matrix Xmat.

Arguments

• Xmat::AbstractArray{T}: state snapshot matrix
• p::Int: order of the Kronecker product

Returns

• kronecker product state snapshot matrix

makeCubicOp(n::Int, inds::AbstractArray{Tuple{Int,Int,Int,Int}}, vals::AbstractArray{Real}, 
which_cubic_term::Union{String,Char}="G") → G or E

Helper function to construct the cubic operator from the indices and values. The indices must be a 1-dimensional array of tuples of the form (i,j,k,l) where i,j,k,l are the indices of the cubic term. For example, for the cubic term $2.5x_1x_2x_3$ for $\dot{x}_4$ would have an index of (1,2,3,4) with a value of 2.5. The which_cubic_term argument specifies which cubic term to construct (the redundant or non-redundant operator). Note that the values must be a 1-dimensional array of the same length as the indices.

Arguments

• n::Int: dimension of the cubic operator
• inds::AbstractArray{Tuple{Int,Int,Int,Int}}: indices of the cubic term
• vals::AbstractArray{Real}: values of the cubic term
• which_cubic_term::Union{String,Char}="G": which cubic term to construct "G" or "E"
• symmetric::Bool=true: whether to construct the symmetric G matrix

Returns

• the cubic operator

makeIdentityCubicOp(n::Int, which_cubic_term::Union{String,Char}="G") → G or E

Helper function to construct the identity cubic operator.

Arguments

• n::Int: dimension of the cubic operator
• which_cubic_term::Union{String,Char}="G": which cubic term to construct "G" or "E"

Returns

• the identity cubic operator

makeQuadOp(n::Int, inds::AbstractArray{Tuple{Int,Int,Int}}, vals::AbstractArray{Real}, 
which_quad_term::Union{String,Char}="H") → H or F or Q

Helper function to construct the quadratic operator from the indices and values. The indices must be a 1-dimensional array of tuples of the form (i,j,k) where i,j,k are the indices of the quadratic term. For example, for the quadratic term $2.5x_1x_2$ for $\dot{x}_3$ would have an index of (1,2,3) with a value of 2.5. The which_quad_term argument specifies which quadratic term to construct. Note that the values must be a 1-dimensional array of the same length as the indices.

Arguments

• n::Int: dimension of the quadratic operator
• inds::AbstractArray{Tuple{Int,Int,Int}}: indices of the quadratic term
• vals::AbstractArray{Real}: values of the quadratic term
• which_quad_term::Union{String,Char}="H": which quadratic term to construct
• symmetric::Bool=true: whether to construct the symmetric H or Q matrix

Returns

• the quadratic operator

make_poly_op(n::Int, inds::AbstractArray{<:NTuple{P,<:Int}}, vals::AbstractArray{<:Real}; 
                nonredundant::Bool=true, symmetric::Bool=true) where P

Helper function to construct the polynomial operator from the indices and values. The indices must be a 1-dimensional array of, e.g., tuples of the form (i,j) where i,j are the indices of the polynomial term. For example, for the polynomial term $2.5x_1x_2$ for $\dot{x}_3$ would have an index of (1,2,3) with a value of 2.5. The nonredundant argument specifies which polynomial operator to construct (the redundant or non-redundant operator). Note that the values must be a 1-dimensional array of the same length as the indices. The symmetric argument specifies whether to construct the operator with symmetric coefficients.

Arguments

• n::Int: dimension of the polynomial operator
• inds::AbstractArray{<:NTuple{P,<:Int}}: indices of the polynomial term
• vals::AbstractArray{<:Real}: values of the polynomial term
• nonredundant::Bool=true: whether to construct the non-redundant operator
• symmetric::Bool=true: whether to construct the symmetric operator

Returns

• the polynomial operator

symmtzrmat(n::Int, p::Int) -> Sp::SparseMatrixCSC{Float64}

Create a symmetrizer matrix of order p for a vector of length n Magnus and Neudecker [1].

Arguments

• n::Int: The length of the vector.
• p::Int: The order of the symmetrizer matrix, e.g., p = 2 for x ⊗ x.

Output

• Sp::SparseMatrixCSC{Float64}: The symmetrizer matrix of order p.

Example

julia> symmtzrmat(2,2)
4×4 SparseMatrixCSC{Float64, Int64} with 6 stored entries:
 1.0   ⋅    ⋅    ⋅ 
  ⋅   0.5  0.5   ⋅
  ⋅   0.5  0.5   ⋅
  ⋅    ⋅    ⋅   1.0

unique_kron_snapshot_matrix(Xmat::AbstractArray{T}, p::Int) where {T<:Number}

Take the p-order unique Kronecker product of each state of the snapshot matrix Xmat.

Arguments

• Xmat::AbstractArray{T}: state snapshot matrix
• p::Int: order of the Kronecker product

Returns

• unique kronecker product state snapshot matrix

unique_kronecker(x::AbstractVector{T}, y::AbstractVector{T}, z::AbstractVector{T}, w::AbstractVector{T}) where T

Unique Kronecker product operation for quadruple Kronecker product.

Arguments

• x::AbstractVector{T}: vector to perform the unique Kronecker product
• y::AbstractVector{T}: vector to perform the unique Kronecker product
• z::AbstractVector{T}: vector to perform the unique Kronecker product
• w::AbstractVector{T}: vector to perform the unique Kronecker product

Returns

• result: unique Kronecker product

Note

This implementation is faster than unique_kronecker_power for p = 4.

unique_kronecker(x::AbstractVector{T}, y::AbstractVector{T}, z::AbstractVector{T}) where T

Unique Kronecker product operation for triple Kronecker product.

Arguments

• x::AbstractVector{T}: vector to perform the unique Kronecker product
• y::AbstractVector{T}: vector to perform the unique Kronecker product
• z::AbstractVector{T}: vector to perform the unique Kronecker product

Returns

• result: unique Kronecker product

Note

This implementation is faster than unique_kronecker_power for p = 3.

unique_kronecker(x::AbstractVector{T}, y::AbstractVector{T}) where T

Unique Kronecker product operation. For example, if

\[x = y = \begin{bmatrix} 1 \\ 2 \end{bmatrix}\]

then

\[\mathrm{unique~kronecker}(x, x) = \begin{bmatrix} 1 \\ 2 \\ 4 \end{bmatrix}\]

Arguments

• x::AbstractVector{T}: vector to perform the unique Kronecker product
• y::AbstractVector{T}: vector to perform the unique Kronecker product

Returns

• result: unique Kronecker product

Note

This implementation is faster than unique_kronecker_power for p = 2.

unique_kronecker_power(x::AbstractArray{T}, p::Int) where T

Unique Kronecker product operation generalized for power p.

Arguments

• x::Union{T,AbstractArray{T}}: vector (or scalar) to perform the unique Kronecker product
• p::Int: power of the unique Kronecker product

Returns

• result: unique Kronecker product

vech(A::AbstractMatrix{T}) → v

Half-vectorization operation. For example half-vectorzation of

\[A = \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix}\]

becomes

\[v = \begin{bmatrix} a_{11} \\ a_{21} \\ a_{22} \end{bmatrix}\]

Arguments

• A: matrix to half-vectorize

Returns

• v: half-vectorized form