Monomials and multiindices

Start from a constructor

Three functions cover the usual shapes. All three return a MultiindexSet in GrLex order.

using MORFE.Multiindices

total_deg = all_multiindices_up_to(2, 4; min_degree = 1) # 1 ≤ |α| ≤ 4
exact_deg = multiindices_with_total_degree(2, 4) # |α| = 4
box = filter(α -> sum(α) ≥ 1, all_multiindices_in_box([3, 3])) # 0 ≤ αᵢ ≤ 3

This code block, considers two variables $z = (z_1, z_2)$ such that a two-component vector of exponents $\alpha = (\alpha_1, \alpha_2)$ represents each monomial $z^\alpha = z_1^{\alpha_1} z_2^{\alpha_2}$.

The first constructor, all_multiindices_up_to, returned the multiindices in 2 variables satisfying $1 \le \alpha_1 + \alpha_2 \le 4$.

The second constructor, multiindices_with_total_degree, returned the multiindices in 2 variables satisfying $\alpha_1 + \alpha_2 = 4$.

The third constructor, all_multiindices_in_box, returned the multiindices in 2 variables satisfying $0 \le \alpha_1, \alpha_2 \le 3$.

See Fig. 01 for a visual representation of the three sets.

Manually define the exponents

The previous constructors are conveniences over one general constructor that takes the exponent vectors directly. It sorts into GrLex and deduplicates for you.

using StaticArrays: SVector

by_hand = MultiindexSet([SVector{2, Int}(a, b)
                        for a in 0:4 for b in 0:4 if 1 ≤ a + b ≤ 4])

by_hand == all_multiindices_up_to(2, 4; min_degree = 1) # true

Two more spellings of the same thing; note the matrix form takes exponents as columns:

MultiindexSet([[1, 0], [0, 1], [1, 1]]) # Vector{Vector{Int}}
MultiindexSet([1 0 1; 
               0 1 1]) # nvars × nmonomials matrix

Combining conditions

In parametric ROM, consider we split the reduced coordinates into master coordinates $z$ and a parameter coordinate $\theta$: we expand to order 4 in $z$ but only to order 2 in $\theta$. Mathematically, we want $\alpha = (\alpha_1, \alpha_2, \alpha_3)$ such that $\alpha_1 + \alpha_2 \le 4$ and $\alpha_3 \le 2$. We also impose the usual $|\alpha| = \alpha_1 + \alpha_2 + \alpha_3 \ge 1$ to remove the constant monomial.

const MAXZ = 4 # max total degree in z
const MAXT = 2 # max degree in θ

# route A — state both conditions in a comprehension
aniso_A = MultiindexSet([SVector{3, Int}(a, b, c)
                        for a in 0:MAXZ for b in 0:MAXZ for c in 0:MAXT
                        if a + b ≤ MAXZ && 1 ≤ a + b + c])

# route B — start from the enclosing box and narrow it
aniso_B = filter(α -> α[1] + α[2] ≤ MAXZ && sum(α) ≥ 1,
                 all_multiindices_in_box([MAXZ, MAXZ, MAXT]))

aniso_A == aniso_B # true — 44 monomials, from a box of 75

Every button below acts on the same lattice: the box $[4,4,2]$ with the origin removed. The red ones are filters and intersect: switch two on and only the monomials satisfying both survive. The green ones add their monomials back on top of whatever the filters left.

Switching on |α| ≤ 4 and α₁ + α₂ ≤ 4 together is the point of this section: the anisotropic cut is not an isotropic expansion of the same order. $|\alpha| \le 4$ in three variables allows $\alpha_3 = 3$ or $4$ while forbidding the total degrees 5 and 6 that the anisotropic set keeps.

Remove multiindices

delete_multiindices returns a new set. There is no in-place variant, so a set already handed to a parametrisation can never be mutated underneath it.

S = all_multiindices_up_to(2, 3; min_degree = 1)

# by explicit exponent — one, a list of them, or another set
delete_multiindices(S, [[3, 0], [0, 3]])
delete_multiindices(S, [1, 1])
delete_multiindices(S, all_multiindices_up_to(2, 1; min_degree = 1))

# by predicate — and `filter` keeps exactly what this drops
odd = delete_multiindices(α -> iseven(sum(α)), S)
even = filter(α -> iseven(sum(α)), S)

length(odd) + length(even) == length(S) # true — they partition S

Superharmonics and spectral radius

Every monomial carries a superharmonic $s(\alpha) = \langle \alpha, \lambda \rangle = \alpha_1 \lambda_1 + \alpha_2 \lambda_2 + \dotsm + \alpha_n \lambda_n$, given the master eigenvalues $\lambda = (\lambda_1, \dotsc, \lambda_n)$ of the linearised system. The superharmonic represents the damping and frequency of that monomial's response. It appears in resonance detection: $\alpha$ is resonant with master mode $r$ when $s(\alpha) \approx \lambda_r $.

We can try to define a spectral cut that keeps only the monomials whose response lands inside a band $|s(\alpha)| < R$ of the spectrum:

λ = [-1 + 2im, -1 - 2im, -1.0 + 0im]
superharmonic(α) = sum(λ .* α)   # s(α) = ⟨λ, α⟩

full = all_multiindices_up_to(3, 4; min_degree = 1)
R = 4.0

band = filter(α -> abs(superharmonic(α)) < R, full)
# 34 monomials → 13

The right-hand panel shows there are only nine distinct points for thirteen monomials: superharmonics collide. $(1,1,0)$ and $(0,0,2)$ both land on $-2$; $(2,1,0)$ and $(1,0,2)$ both land on $-3+2i$ by different routes. Collisions are the reason resonance detection exists.

This set is not downward closed: it keeps $\alpha = (2,1,0)$ while dropping its factor $(2,0,0)$; and keeps $(1,2,0)$ while dropping $(0,2,0)$. Step 06 explains why this is a problem.

Downward closure

To compute a monomial $z^\alpha$, we require all its factors to be computed first. The factors of $\alpha$ are all $\beta$ in the hyperrectangle (the box) satifying $\beta \le \alpha$ component-wise, i.e. $\beta_i \le \alpha_i$. Sets safisfying this property are called downward closed.

E.g. $z_1^2 z_2^1 z_3^0$ has factors $\{z_1, z_2, $ $z_1^2$ $, z_1 z_2, z_1^2 z_2\}$. Therefore $\alpha=(2,1,0)$ has factors $\{(1,0,0), (0,1,0),$ $(2,0,0)$ $, (1,1,0), (2,1,0)\}$.

And this explains why the MultiindexSet from Step 05 is not downward closed:

is_downward_closed(band) # false
# [2,1,0] is kept (|s| = 3.606), but not its factor [2,0,0] with |s| = 4.472 ≥ R
# [1,2,0] is kept (|s| = 3.606), but not its factor [0,2,0] with |s| = 4.472 ≥ R

The fix is to take the downward closure — add back every factor of every member. The result is legal, and still far smaller than the full expansion:

function downward_closure(S::MultiindexSet{N}) where {N}
    closed = Set{SVector{N, Int}}()
    for α in S.exponents
        for idx in CartesianIndices(ntuple(i -> 0:α[i], N))
            push!(closed, SVector{N, Int}(Tuple(idx)))
        end
    end
    delete!(closed, zero(SVector{N, Int})) # DPIM needs min_degree ≥ 1
    return MultiindexSet(collect(closed))
end

is_downward_closed(downward_closure(band)) # true — 15 monomials, vs 34
# the closure added exactly [2,0,0] and [0,2,0]

The full contract for a custom mset:

All five are enforced by validate_multiindex_set, which parametrise and solve_parametrisation both run before solving. You can check a set the moment you build it, and it names the offending exponent:

validate_multiindex_set(band, 3, 3)
# ArgumentError: custom mset is not downward closed: [2, 1, 0] is a member
# but its factor [2, 0, 0] is not. The graded solve reads W[α - β + eᵢ] and
# factorises α over members of mset, so a missing factor is read as zero
# and silently corrupts the right-hand side.

Pass validate_mset = false to skip the check if you have already done it yourself. Then hand it over to parametrise to get the reduced order model:

W, R = parametrise(model, spectral_data, mset)