SOAR, step by step

A frequency-domain finite element model with frequency-independent matrices, time convention \(e^{\jj\omega t}\), \(\omega=2\pi f\):

\(\vu\in\mathbb C^N\) is the nodal pressure, \(\vF\) the load, and \(H\) a linear output (one row \(\vL^{\Tr}\), no conjugation).

The wall impedance is a dial between a smooth response and sharp resonances (rigid-tube resonances at \(f_n=nc/2L=171.5\,n\) Hz):

Each frequency of the direct sweep assembles \(\vZ(\omega)\), factorises it and solves. The factorisation has to be redone at every frequency, because \(\vZ\) changes with \(\omega\). Model order reduction replaces \(\vu(\omega)\) by \(\vV\vy(\omega)\) with one frequency-independent basis \(\vV\in\mathbb C^{N\times n}\), \(n\ll N\). Moment matching chooses \(\vV\) from the Taylor expansion of \(\vu\) at one frequency \(\omega_0\). SOAR is a stable way to compute that basis when \(\vZ\) is quadratic.

Four properties of (1.1) are used throughout:

• Quadratic in \(\omega\) with constant matrices. A frequency-dependent impedance \(Z_s(\omega)\) breaks this (Section 12).
• Complex symmetric. \(\vK,\vC,\vM\) are real symmetric, so \(\vZ(\omega)^{\Tr}=\vZ(\omega)\), but \(\vZ^{\Hm}\neq\vZ\).
• \(\vK\vone=\vzero\). No Dirichlet boundary, so \(\vZ(0)=\vK\) is singular: \(\omega=0\) is a pole and \(\omega_0=0\) is not allowed. Physically, an acceleration source has an unbounded static displacement.
• \(\vC\) has one entry. The damping is not proportional (\(\vC\neq\alpha\vK+\beta\vM\)), so the real modes of \((\vK,\vM)\) do not decouple the system. This is the setting where moment matching is the natural tool.

Put \(\omega=\omega_0+\delta\) and expand:

• The split is exact. \(\vZ\) is a polynomial of degree 2, so \(\vZ_1=\vZ'(\omega_0)\), \(\vZ_2=\tfrac12\vZ''(\omega_0)\) and no term is dropped.
• Only \(\vZ_0\) is factorised. \(\vZ_1\) and \(\vZ_2\) are only multiplied with vectors.
• \(\vZ_0\) must be invertible, so \(\omega_0\) must not be a pole. All poles of the tube satisfy \(\mathrm{Im}\,\omega\ge0\) (Section 8), and the only real one is \(\omega=0\) (Section 3). Any real \(\omega_0>0\) therefore works in exact arithmetic. Near a lightly damped resonance \(\vZ_0\) is close to singular, but sparse LU with pivoting stays backward stable.

# soar_tube/moments.py and soar_tube/fom.py
Z0, Z1, Z2 = expansion_matrices(K, C, M, omega0)   # K + 1j*omega0*C - omega0**2*M,  1j*C - 2*omega0*M,  -M
solve0 = lu_solver(Z0)                            # one sparse LU (SciPy SuperLU), returns solve(b)

Because \(\vZ_0\) is invertible, \(\vu(\omega)=\vZ(\omega)^{-1}\vF\) is analytic in a disk around \(\omega_0\). Its radius \(\rho\) is the distance from \(\omega_0\) to the nearest pole, i.e. the nearest root of \(\det\vZ(\omega)=0\):

Insert (2.1) and (3.1) into \(\vZ\vu=\vF\) and collect equal powers of \(\delta\):

Hence the three-term recurrence that defines everything that follows:

Each new moment costs two sparse products and one forward and backward substitution with the LU factors of \(\vZ_0\). \(\vA\) and \(\vB\) are dense and are never formed. The output moments are \(m_i=\vL^{\Tr}\vr_i\), so \(H(\omega_0+\delta)=\sum_i m_i\delta^i\).

The tube's poles are known in closed form. With \(R=(Z_s-\rho c)/(Z_s+\rho c)\), \(\det\vZ=0\) for the continuous problem when \(R\,e^{-2\jj kL}=1\):

all with the same damping (analytic_poles_hz). The FE poles agree to …. For the moderate wall, the pole nearest to \(f_0=1040\) Hz is at …. So (3.1) converges only within … of 1040 Hz, i.e. on …. Yet the SOAR model with 20 vectors is accurate to 0.1 % over 20–2000 Hz. A reduced model is not a truncated Taylor series. It is a rational function whose poles move towards the true poles (Section 8). This is the main conceptual point of moment matching.

Asymptotic waveform evaluation (AWE, Pillage & Rohrer 1990) builds a Padé approximant from the scalar moments \(m_i\). One could also project on \(\spn\{\vr_0,\dots,\vr_{n-1}\}\) computed with (3.2). Both fail after a handful of moments in double precision (… on the tube). The reason is that (3.2) is a power iteration.

Let \(\vZ(\omega_p)\vx=\vzero\), \(\omega_p=\omega_0+\delta_p\). Then \((\vZ_0+\delta_p\vZ_1+\delta_p^2\vZ_2)\vx=\vzero\). Multiply by \(-\vZ_0^{-1}/\delta_p^2\) and set \(\mu=1/\delta_p\):

So every pole \(\omega_p\) gives an eigenvalue \(\mu=1/(\omega_p-\omega_0)\) of \(\vH\). Written for pairs of moments, (3.2) is

This is the power method on \(\vH\). The normalised moments turn towards the eigenvectors with the largest \(|\mu|\), i.e. the poles closest to \(\omega_0\), and each new moment adds less and less new direction. Scaling does not help: it cures the growth of \(\|\vr_i\|\sim\rho^{-i}\), not the loss of direction.

Definition (Bai & Su 2005). For \(\vA,\vB\in\mathbb C^{N\times N}\) and \(\vr_0\neq\vzero\):

With \(\vA,\vB,\vr_0\) from (3.2), \(\cG_n\) is exactly the span of the first \(n\) moments of \(\vu\). For \(\vB=\vzero\) it reduces to the ordinary Krylov subspace \(\cK_n(\vA;\vr_0)\).

Linearisation. Stack two consecutive moments, \(\vw_i=\begin{bmatrix}\vr_i\\\vr_{i-1}\end{bmatrix}\in\mathbb C^{2N}\). The top block of \(\vH\vw_{i-1}\) is \(\vA\vr_{i-1}+\vB\vr_{i-2}=\vr_i\), and the bottom block is \(\vr_{i-1}\). So

The second-order Krylov subspace is the top half of an ordinary Krylov subspace of dimension \(2N\). By (4.1), the eigenvalues of \(\vH\) are \(\mu=1/(\omega_p-\omega_0)\). \(\vH\) is the shift-and-invert operator of a companion linearisation, which is why Krylov methods on \(\vH\) find the poles nearest \(\omega_0\) first.

This gives two routes:

• Route A, linearise then reduce. Arnoldi on \(\vH\) gives an orthonormal \(\vV_n\in\mathbb C^{2N\times n}\). The vectors are twice as long, so memory and Gram–Schmidt cost double. The top block spans \(\cG_n\) but is not orthonormal and must be re-orthonormalised. Projecting the first-order (linearised) system instead gives a first-order ROM that loses the \(\vK,\vC,\vM\) structure and the stability argument of Section 8.
• Route B, SOAR. Orthonormalise only the top block. Carry the bottom block as auxiliary vectors \(\vp_j\) of length \(N\). Never form a vector of length \(2N\).

Goal. Build \(\vQ_n=[\vq_1,\dots,\vq_n]\) with \(\vQ_n^{\Hm}\vQ_n=\vI\) and \(\spn\{\vq_1,\dots,\vq_j\}=\cG_j\) for every \(j\le n\), using vectors of length \(N\) only.

Ansatz (Bai & Su). Look for auxiliary vectors \(\vP_n=[\vp_1,\dots,\vp_n]\) and an \((n+1)\times n\) upper Hessenberg matrix \(\That_n=(t_{ij})\) such that the stacked vectors satisfy an Arnoldi-type relation with \(\vH\):

The difference with Arnoldi is that only the top block \(\vQ\) is required to be orthonormal. Write (6.1) block row by block row:

Column \(j\) of the top equation. \(\vA\vq_j+\vB\vp_j=\sum_{i=1}^{j+1}\vq_i\,t_{ij}\). Multiply by \(\vq_k^{\Hm}\), \(k\le j\), and use orthonormality:

This is Gram–Schmidt of the candidate \(\vA\vq_j+\vB\vp_j\) against \(\vQ_j\).

Column \(j\) of the bottom equation. \(\vq_j=\sum_{i=1}^{j+1}\vp_i\,t_{ij}\), so

The bottom block is not orthogonalised. It receives the same coefficients \(t_{ij}\), because (6.1) requires it. Start: \(\vw_0=[\vr_0;\vzero]\) scaled by \(1/\|\vr_0\|\), so \(\vq_1=\vr_0/\|\vr_0\|\) and \(\vp_1=\vzero\).

Meaning of \(\vp_j\). The recurrence \(\vr_i=\vA\vr_{i-1}+\vB\vr_{i-2}\) needs two vectors, the current moment and the previous one. SOAR does the same with \(\vq_j\) and its companion \(\vp_j\), which stands for "the previous vector" in the stacked pair \([\vr_i;\vr_{i-1}]\).

Taking the first \(n\) rows of \(\That_n\) as \(\vT_n\), (6.1) reads

the analogue of the Arnoldi decomposition. Bai & Su introduced SOAR for the quadratic eigenvalue problem. For eigenvalues they project the QEP onto \(\vQ_n\) (Rayleigh–Ritz) rather than using \(\vT_n\). The reduced model of Section 8 is that same projection.

If \(t_{j+1,j}=0\) in (6.4), the candidate brings no new direction to the top block and the recurrence cannot divide by it. Bai & Su's Algorithm 2 separates two cases:

• Deflation. \(\vA\vq_j+\vB\vp_j\in\spn\vQ_j\), but the stacked vector \([\vzero;\vs]\) is still new. Set \(t_{j+1,j}=1\), \(\vq_{j+1}=\vzero\), \(\vp_{j+1}=\vs\), and continue. \(\cG\) pauses while \(\cK(\vH;\vw_0)\) keeps growing, and later steps can bring new directions back into the top block. The zero columns of \(\vQ\) are dropped at the end.
• Breakdown. If also \(\vs\in\spn\{\vp_i:\ \vq_i=\vzero\}\), then \(\cK_j(\vH;\vw_0)\) is invariant under \(\vH\). Every moment lies in \(\cG_j\), so \(\vu(\omega)\in\cG_j\) on the Taylor disk, and by analytic continuation for all \(\omega\). In exact arithmetic the ROM is then exact (a "lucky" breakdown) and the procedure stops.

Approximate \(\vu(\omega)\approx\vV\vy(\omega)\) and ask the residual to be orthogonal to \(\spn\vV\) (Galerkin): \(\vV^{\Hm}\big(\vF-\vZ(\omega)\vV\vy\big)=\vzero\). This gives

and \(H_r(\omega)=\vL_r^{\Tr}\vy(\omega)\) with \(\vL_r^{\Tr}=\vL^{\Tr}\vV\) (project, rom_solve, rom_sweep).

• Same structure. The ROM is quadratic in \(\omega\) with constant \(n\times n\) matrices: \(\mathcal O(n^3)\) per frequency and no reassembly.
• Nested. \(\vV_k\) is the first \(k\) columns of \(\vV_n\). The size-\(k\) ROM is therefore the leading \(k\times k\) block of (8.1): one projection serves every size (truncate).
• Only the subspace matters. Replacing \(\vV\) by \(\vV\vS\) (\(\vS\) invertible) changes \(\vy\) to \(\vS^{-1}\vy\) and leaves \(H_r\) unchanged. In exact arithmetic, SOAR, Arnoldi on \(\vH\) with a QR of the top block, the explicit-moment basis and WCAWE (for this quadratic \(\vZ\)) all give the same ROM. They differ only in floating point.
• The basis depends on \(\vF\), not on \(\vL\). One basis serves every output: the wall pressure and the space-averaged pressure need the same \(n\) (Section 10).

The proof only uses \(\cG_n\subseteq\spn\vV\). The job of SOAR is to deliver that subspace with an orthonormal, well-conditioned basis.

What kind of approximant is \(H_r\)?

\(\det\vZ_r(\omega)\) is a polynomial of degree \(2n\), so \(H_r\) is rational with \(2n\) poles and a numerator of degree at most \(2n-2\). A Padé approximant with that many poles would match about \(4n\) moments. The one-sided ROM matches \(n\), the two-sided one \(2n\), and the rest of the freedom goes to keeping \(\vK,\vC,\vM\) and stability. The ROM poles lock onto the true poles nearest \(\omega_0\) first, so the accurate band grows outward from \(f_0\) as \(n\) increases. The live demo shows this.

Orthogonality

Classical Gram–Schmidt applied twice (CGS2, "twice is enough") keeps \(\vV\) orthonormal to working precision: \(\|\vV^{\Hm}\vV-\vI\|=\)… at 1040 Hz (moderate wall), and the recurrence (6.1) holds to … (check_basis, soar_relation_residual). A single CGS pass loses orthogonality as the basis converges.

The auxiliary vectors grow

By (6.7), \(\|\vp_{j+1}\|=\|\Ts^{-1}\ve_j\|\). The diagonal of \(\Ts\) holds the subdiagonal entries \(t_{k+1,k}\), the size of each new direction. As the subspace captures the poles near \(\omega_0\), these entries shrink relative to the other \(t_{ik}\). This is the same mechanism that makes the Arnoldi subdiagonal shrink when Ritz values converge. \(\Ts^{-1}\) then grows geometrically.

The consequence: the candidate \(\vA\vq_j+\vB\vp_j\) is dominated by \(\vB\vp_j\), and that part lies almost entirely in \(\spn\vQ_j\). Gram–Schmidt removes it. What remains is the fraction

of the candidate. A rounding error of order \(\varepsilon\|\vA\vq_j+\vB\vp_j\|\) in the candidate becomes a relative error of about \(\varepsilon/\eta_j\) in \(\vq_{j+1}\). When \(\eta_j\) reaches the tolerance \(10^{-12}\), line 7 of the algorithm declares a deflation: the new direction is indistinguishable from round-off.

Numerical deflations, then SOAR stops

On the tube, every run of 40 vectors ends the same way. Once \(\|\vp_j\|\) reaches \(10^{14}\)–\(10^{16}\), \(\eta_j\) falls below the tolerance, line 8 deflates, and the following candidates are dominated by round-off as well. Left alone, the deflations would come in a cascade until the span test declares a breakdown, without adding a single vector; soar_basis therefore stops after three deflations in a row (line 9, stopped_reason = "deflations"). The basis then stops growing. These are numerical events, not the exact deflation of Section 7, and they decide how far SOAR can go:

• At a good expansion point it does no harm. At 1040 Hz the basis stops at … vectors, when the error is already ….
• Far from part of the band it does. At 200 Hz, at the edge of a 20–2000 Hz band, the basis stops at … vectors with a … error. A stable basis of the same subspace (Arnoldi on the \(2N\) linearisation, then QR of the top block) keeps converging, below. The stall is a precision effect, not a property of the subspace.

Rescaling \(\omega\) does not help. With \(\delta=s\delta'\): \(\vA\to s\vA\), \(\vB\to s^2\vB\), every \(\vp_j\to\vp_j/s\) and every \(t_{ij}\to s\,t_{ij}\). \(\eta_j\) does not change. The cause is structural: SOAR normalises the top block only, so it divides by \(t_{j+1,j}\) even when the stacked vector \([\vq;\vp]\) still carries plenty of new information. TOAR (Section 12) normalises the whole stacked vector and avoids the division.

In practice: always use CGS2, monitor \(\eta_j\), \(\|\vp_j\|\) and the deflations reported by soar_basis, and when they appear before the error is small enough, either switch to TOAR or add a second expansion point instead of pushing \(n\).

How many vectors, without the full model

• Nested difference. All sizes \(k\le n\) come for free (leading blocks), so \(\|H_n-H_{n-d}\|/\|H_n\|\) over the band (\(d=2\) to 5) is a free error indicator. Stop when it falls below the tolerance.
• Residual at \(\mathcal O(n^2)\) per frequency. \(\vR(\omega)=\vF-\vZ(\omega)\vV\vy=\vW\vc(\omega)\) with \(\vW=[\vF,\ \vK\vV,\ \vC\vV,\ \vM\vV]\) and \(\vc=[1;\,-\vy;\,-\jj\omega\vy;\,\omega^2\vy]\). Precompute \(\vG=\vW^{\Hm}\vW\) once (size \(3n+1\)); then \(\|\vR(\omega)\|^2=\vc^{\Hm}\vG\vc\) costs nothing that depends on \(N\). Since error \(\le\|\vZ(\omega)^{-1}\|\,\|\vR(\omega)\|\), a small residual can still hide a large error at a lightly damped peak.
• Rule of thumb: about two vectors per mode in the band (tube: 11 resonances, \(n\approx19\)–20 for \(10^{-3}\); car of Section 13: 8 modes, \(n\approx20\)).

One point or several

Multi-point SOAR concatenates the bases built at several \(\omega_k\), orthonormalises them and projects once (multipoint_basis, union_basis). Theorem 8.1 applies at every \(\omega_k\), because each \(\cG_{n_k}(\omega_k)\) lies in the union. Each point costs one LU. Counted in units of one SOAR vector, with \(\beta\) = cost of one LU / cost of one vector:

Extra points save one or two vectors on the tube, because one point already covers the band. With \(\beta\approx2\) in 1D that is a wash; with \(\beta\approx72\) on a 3D model it is a clear loss. An extra point pays only if it saves more vectors than one LU costs, or if the factorisations run in parallel.

Where the time goes

… In 1D a sparse LU costs about as much as a back-substitution, so the speed-up is modest. When the LU dominates, as in 3D, the speed-up approaches the number of frequencies: the whole SOAR model costs about as much as one frequency of the direct sweep (Section 13).

TOAR in a few lines

The bottom block of \(\vH\vv\) is the top block of \(\vv\). So every block of every Arnoldi vector of \(\cK_j(\vH;\vw_0)\) lies in \(\cG_{j+1}=\spn\{\vr_0,\dots,\vr_j\}\). Write the orthonormal Arnoldi basis as

One step computes \(\vr=\vA(\vU\vs_1)+\vB(\vU\vs_2)\) (one substitution), orthogonalises \(\vr\) against \(\vU\) (first level, length \(N\)), and then orthogonalises the new coordinate vector against \(\vS\) (second level, small). Memory is about \(N(n+1)\), like SOAR. The stacked basis is orthonormal, like Arnoldi on \(\vH\), so nothing is divided by a small top-block norm. The subspace and the ROM are the same as SOAR's. TOAR is available in SLEPc (PEPTOAR, Campos & Roman 2016) and extends to higher-degree polynomials (Kressner & Roman 2014).

Other directions

• Frequency-dependent load. For \(\vF(\omega_0+\delta)=\sum_k\vF_k\delta^k\), the recurrence becomes \(\vZ_0\vr_i=\vF_i-\vZ_1\vr_{i-1}-\vZ_2\vr_{i-2}\). For a velocity source \(\jj\omega\rho\,v\), it is simpler to solve with the constant part and multiply the output by \(\jj\omega\) (the problem is linear). This is why the tube is driven by an acceleration.
• Frequency-dependent impedance or admittance. \(\vZ\) is no longer quadratic. Use WCAWE (Souza Lenzi et al. 2013, §6.1.2), or fit a rational impedance and add auxiliary unknowns to recover a polynomial system.
• Several loads. Block SOAR with \(\vF=[\vF_1,\dots,\vF_m]\).
• Two-sided projection. \(2n\) moments with a second, output-dependent subspace (Theorem 8.2 is the symmetric special case).
• Eigenvalues. SOAR was designed for the QEP. Rayleigh–Ritz on \(\vV\) gives the ROM poles shown in the demo.

The repository provides the car case without any FE code:

• data/car_cavity/car_cavity.md: geometry, boundary conditions, parameters and what was assumed.
• car_cavity.step: the cavity solid (any CAD tool or gmsh); car_cavity_geometry.py: the gmsh script that builds and labels it.
• car_cavity_p2_19072dofs.npz: assembled \(\vK,\vC,\vM,\vF\) and the ear outputs of a converged 19 072-dof P2 model, ready for soar_basis; examples/05_car_cavity.py runs SOAR on it.