Pencil Decomposition: Methods & Applications

• Pencil decomposition is a versatile concept encompassing parallel grid layouts, spectral methods, tensor rank decompositions, and structural reductions of matrix/operator pairs.
• In high-performance computing, it enables efficient data distribution by splitting multidimensional grids into pencils, reducing memory footprint and communication overhead for large-scale simulations.
• Applications extend to dynamic lung MRI and optimal control, where pencil-based methods separate signal components and isolate system constraints to enhance analysis and reproducibility.

In the cited literature, pencil decomposition does not denote a single method. It refers, depending on context, to a two-dimensional domain decomposition for large parallel grids, to a matrix pencil decomposition for spectral analysis of sampled signals, or to a reduction or block decomposition of matrix/operator pencils used to expose spectrum, constraints, or canonical structure [(Kawazura, 2022); (Pusterla et al., 2024); (Verdier, 2012)]. The term is therefore field-dependent: in large-scale simulation it is primarily a data-layout strategy; in signal processing it is a generalized-eigenvalue-based spectral method; and in linear systems, differential-algebraic equations, and operator theory it is a structural decomposition of pairs such as $(E,A)$ or $\lambda B-A$.

1. Pencil decomposition as a two-dimensional parallel grid layout

In high-performance PDE solvers, pencil decomposition is a parallel data layout in which a three-dimensional grid is split across MPI tasks in two directions at once, leaving one direction local to each process. The standard contrast is with slab decomposition, where the grid is distributed in only one direction. For an $N^3$ grid, slab decomposition allows at most $N$ MPI processes, whereas pencil decomposition allows up to $N^2$ MPI processes in principle (Kawazura, 2022).

This distinction is decisive for methods dominated by multidimensional transforms. In the pseudospectral method, the equations are evolved in Fourier space, while physical space is used only to compute nonlinear products. Computation therefore alternates between inverse FFTs, nonlinear term evaluation, and forward FFTs. For large three-dimensional domains, the FFTs dominate runtime, and process-to-process transposes are often more expensive than the FFT arithmetic itself. Pencil decomposition reduces the per-process memory footprint and increases concurrency, which is important for strong scaling in turbulence simulations that require high resolution to resolve the inertial range (Kawazura, 2022).

A closely related use appears in a finite-difference DNS code for wall-bounded turbulence. There, each MPI process owns a thin three-dimensional subdomain shaped like a pencil, with one direction complete on a process and the other two split among process rows and columns. For the wall-bounded flows considered, the preferred orientation aligns the pencil with the wall-normal direction $x$, so that each process owns the full $x$-extent for a subset of $(y,z)$ planes. This orientation is essential because the wall-normal viscous terms remain implicit, and the associated tridiagonal solves can then be carried out locally for each $(y,z)$ pair without all-to-all communication (Poel et al., 2015).

The numerical motivation in the wall-bounded case is anisotropy of both geometry and timestep constraints. The code keeps the wall-normal viscous terms implicit, makes the wall-parallel viscous terms explicit, and relies on the fact that at high driving the CFL condition is more restrictive than the explicit viscous constraint in the wall-parallel directions. A direct consequence is that implicit integration of wall-parallel viscous terms is no longer required, avoiding non-local information exchange for those derivatives (Poel et al., 2015).

2. FFT workflows, transpositions, and scalability

The FFT-oriented implementation of pencil decomposition is illustrated by the pseudospectral MHD code Calliope. It adopts the P3DFFT library for parallel 3D FFTs using pencil decomposition, uses a stride-1 data structure, and stores local arrays in a layout that is localized in one dimension, with FFTs and transposes managed by P3DFFT (Kawazura, 2022). The local array layout is given as

$$\begin{array}{c} \text{Physical space}\ (N^{(l)}_y,\, N^{(l)}_z/M_1,\, N^{(l)}_x/M_2) \end{array} \;\; \Longleftrightarrow \;\; \begin{array}{c} \text{Fourier space}\ (N^{(k)}_x,\, N^{(k)}_z/M_2,\, N^{(k)}_y/M_1) \end{array}$$

with

$\lambda B-A$0

and $\lambda B-A$1 the total number of MPI processes. In this layout, the $\lambda B-A$2-direction is localized in physical space and the $\lambda B-A$3-direction is localized in Fourier space (Kawazura, 2022).

That choice is tied to the shearing-box remapping procedure used for local magnetorotational turbulence. In shearing coordinates, the radial Fourier wavenumber evolves as

$\lambda B-A$4

and remapping is performed periodically at

$\lambda B-A$5

The fields are rearranged as

$\lambda B-A$6

after which $\lambda B-A$7 is reset to $\lambda B-A$8. Because $\lambda B-A$9 is localized in Fourier space, this rearrangement can be done without transferring data between processors during remapping (Kawazura, 2022). The paper states that Calliope is currently the only pseudospectral code that computes magnetorotational turbulence using pencil decomposition (Kawazura, 2022).

The wall-bounded finite-difference code exhibits a different but related transpose pattern. The pressure correction is solved through 2D FFTs in the periodic directions and a tridiagonal solve in the wall-normal direction. The algorithm transposes data through the sequence $N^3$0, and then reverses that sequence after the tridiagonal solve. The non-neighbor transposes, especially $N^3$1, are implemented using MPI_ALLTOALLW, while the others use MPI_ALLTOALLV (Poel et al., 2015).

Step	Slab A2A	Pencil A2A
$N^3$2 computation	6	0
Pressure correction	2	6
Scalar equation	2	0
Total	10	6

This communication count summarizes the effect of the semi-implicit pencil-decomposed formulation: all-to-all communication is eliminated from the $N^3$3 and scalar updates, leaving only the six transposes required by the pressure Poisson solve (Poel et al., 2015).

The reported scaling results make the architectural consequences explicit. For incompressible MHD on a $N^3$4 grid, Calliope maintains almost ideal scaling up to $N^3$5 cores on Fugaku (Kawazura, 2022). The finite-difference wall-bounded code shows near-linear strong scaling up to 32K cores for $N^3$6, slight loss at 64K cores, and weak scaling close to linear from 2 up to 16K cores with $N^3$7 points per core (Poel et al., 2015).

3. Matrix pencil decomposition in dynamic lung MRI

In functional lung MRI, matrix pencil decomposition is a spectral-analysis method rather than a domain decomposition. The TrueLung pipeline applies it to free-breathing, contrast-agent-free pulmonary MRI based on 2D+t coronal ufSSFP image series. The pipeline includes data transfer and sorting, initial quality checks, respiratory phase selection, image registration, neural-network whole-lung and lobar segmentation, voxel-wise matrix pencil functional analysis, defect quantification, and reporting (Pusterla et al., 2024).

The central computational role of matrix pencil decomposition in this setting is to separate the observed signal-time series into components associated with respiration, cardiac perfusion, and residual or unwanted components. The paper presents it as an improved variant of Fourier decomposition that is more robust than Fourier-based methods to irregular breathing, variable cardiac frequency, and time-series truncation (Pusterla et al., 2024). This robustness is material in free-breathing acquisitions, where strict periodicity is absent.

The signal model is the standard damped-exponential form

$N^3$8

where $N^3$9 are complex amplitudes and $N$0. Classical matrix-pencil matrices are then built from the sampled time series as Hankel-like data matrices $N$1 and $N$2, and the poles $N$3 are estimated through a generalized eigenvalue problem or an equivalent matrix pencil decomposition (Pusterla et al., 2024). In the TrueLung workflow, the signal intensity in the segmented lung is first integrated, spectral analysis using matrix pencil detects global signal modulations, these modulations are classified as respiratory or cardiac using k-means clustering and reference frequency ranges, and the corresponding amplitudes and phases are retrieved by matrix-pencil decomposition and least-squares fitting (Pusterla et al., 2024).

The outputs are fractional ventilation, perfusion, and blood arrival time maps. These maps are then thresholded to derive the defect percentages RFv and RQ. The study used a threshold equal to 75% of the median value from each voxel distribution, while the appendix notes the general practice of using 60–80% of the median depending on the application. For ventilation quantification, vessels are removed because they behave like non-ventilated regions in ventilation maps; the paper identifies vessels from perfusion maps using the 85th percentile of signal intensities (Pusterla et al., 2024).

The reported validation is strongly oriented toward pipeline reproducibility and the effect of automated segmentation. Functional imaging was performed at 7.9 $N$4 1.8 coronal slice positions per patient, with a total scan time of 6min 20s per patient on average. The whole pipeline required 20min calculation time per subject, with registration taking 17 min 6 s, segmentation 2 min 10 s, matrix pencil analysis 30 s, and defect quantification plus report generation < 1 s (Pusterla et al., 2024). Automated whole-lung quantification was satisfactory in 88% of patients and 97% of slices, while lobar quantification was satisfactory in 73% of patients and 93% of slices. The relative differences in RFv and RQ between fully automated and manually refined data were reported as marginal (Pusterla et al., 2024).

4. Pencil-based tensor rank decomposition

In tensor analysis, a pencil decomposition can denote a tensor rank decomposition algorithm that first reduces a third-order tensor to an $N$5 tensor, converts that compressed tensor into a linear matrix pencil, and then recovers factor matrices through a generalized eigendecomposition (Beltrán et al., 2018). The setting is the canonical polyadic decomposition (CPD)

$N$6

A pencil-based algorithm first projects along the third mode using a matrix $N$7 with orthonormal columns and forms

$N$8

After further orthogonal Tucker compression, the paper obtains a core tensor $N$9 with slices

$N^2$0

If $N^2$1 is nonsingular, then

$N^2$2

so the pair $N^2$3 is the matrix pencil and the eigenvectors recover $N^2$4. The generalized eigenvalue problem is

$N^2$5

The paper’s main result is negative: for $N^2$6, every pencil-based algorithm is numerically forward unstable on some open set of rank-$N^2$7 tensors (Beltrán et al., 2018). The stated cause is not unstable linear-algebra subroutines, but a conditioning mismatch introduced by the projection to two slices. The condition number of the CPD can be much larger for $N^2$8 tensors than for the original $N^2$9 tensor. The paper highlights that the projected vectors $x$0 may become nearly aligned even when the original $x$1 are not, and gives the lower bound

$x$2

The instability is therefore structural. The paper proves that the excess factor $x$3 can be made arbitrarily large, establishes lower bounds for the limiting distribution of the condition number of random CPDs, and states that for random tensor rank decompositions one should anticipate a loss of precision of a few digits (Beltrán et al., 2018). It also recommends Newton-based refinement when high accuracy is required.

5. Reduction theory for matrix and operator pencils

A separate and extensive line of work uses decomposition to analyze matrix pencils and operator pencils as pairs of maps rather than as sampled-signal models. In one formulation, the pencil is a pair of bounded linear operators

$x$4

on Banach spaces, with reduction procedures designed to remove constraints and expose the less implicit part of equations such as $x$5 or $x$6 (Verdier, 2012). The two basic operations are observation reduction and control reduction.

Observation reduction is defined through

$x$7

together with reduced operators $x$8 and the quotient map

$x$9

which is the pivot operator for the observation reduction. Control reduction uses

$x$0

with induced maps $x$1 and pivot operator

$x$2

These reductions are conjugate to one another, and under the paper’s normality assumptions they commute canonically:

$x$3

with

$x$4

The same framework gives an intrinsic criterion for irreducibility:

$x$5

It also characterizes the resolvent set

$x$6

showing that invertible pivot operators preserve the resolvent under reduction, while noninvertible pivot operators force the resolvent to be empty (Verdier, 2012).

In finite-dimensional vector spaces, an invariant reduction theory for matrix pencils

$x$7

recovers the classical Kronecker and Weierstraß pictures without choosing bases at the outset (Verdier, 2012). The intrinsic reduction is

$x$8

iterated as $x$9. The associated quotient maps define the defect sequences

$(y,z)$0

interpreted respectively as constraint, observation, and control defects. The paper states that a pencil is regular iff

$(y,z)$1

and that the reduction framework recovers the Kronecker blocks $(y,z)$2, $(y,z)$3, and their duals from these defect sequences (Verdier, 2012).

This reduction picture extends beyond vector spaces. In an abelian category, one considers a system

$(y,z)$4

with auxiliary morphisms

$(y,z)$5

Observation and control reductions become

$(y,z)$6

and the reductions still commute (Verdier, 2016). The resulting invariant sequences split into nilpotent, observation, and control defects; in the vector-space case these defect objects become the multiplicities of the classical Kronecker blocks (Verdier, 2016).

6. Specialized pencil decompositions in spectral analysis and control

Several applications use pencil decomposition to obtain tractable block structures for spectral or control-theoretic analysis. One example is the one-dimensional heat equation with mixed boundary conditions, rewritten as a non-self-adjoint $(y,z)$7 unbounded block operator matrix pencil

$(y,z)$8

After setting $(y,z)$9 and introducing

$(y,z)$0

the problem becomes

$(y,z)$1

on

$(y,z)$2

The pencil is non-self-adjoint because of the derivative coupling in the block structure, so pseudospectra rather than eigenvalues alone are used to assess sensitivity (G. et al., 2023). The generalized eigenvalues are

$(y,z)$3

and the paper develops spectral and pseudospectral enclosures through generalized Frobenius–Schur factorizations of $(y,z)$4 block operator pencils (G. et al., 2023). For the heat problem, the specialization

$(y,z)$5

shows how the non-self-adjoint pencil is controlled by a self-adjoint reference operator plus a coupling term (G. et al., 2023).

Another specialized use appears in the finite-horizon LQ problem with two-sided boundary conditions. There, the first-order optimality system yields the extended symplectic pencil

$(y,z)$6

with

$(y,z)$7

If $(y,z)$8 solves the constrained generalized discrete-time algebraic Riccati equation, the paper proves the existence of invertible matrices $(y,z)$9 such that

$$\begin{array}{c} \text{Physical space}\ (N^{(l)}_y,\, N^{(l)}_z/M_1,\, N^{(l)}_x/M_2) \end{array} \;\; \Longleftrightarrow \;\; \begin{array}{c} \text{Fourier space}\ (N^{(k)}_x,\, N^{(k)}_z/M_2,\, N^{(k)}_y/M_1) \end{array}$$0

This decomposition does not require regularity of the extended symplectic pencil, modulus controllability, or solvability of a closed-loop Lyapunov equation (Ferrante et al., 2012). In the singular case, the finite generalized eigenvalues come from the dynamics of $$\begin{array}{c} \text{Physical space}\ (N^{(l)}_y,\, N^{(l)}_z/M_1,\, N^{(l)}_x/M_2) \end{array} \;\; \Longleftrightarrow \;\; \begin{array}{c} \text{Fourier space}\ (N^{(k)}_x,\, N^{(k)}_z/M_2,\, N^{(k)}_y/M_1) \end{array}$$1 induced on the quotient space $$\begin{array}{c} \text{Physical space}\ (N^{(l)}_y,\, N^{(l)}_z/M_1,\, N^{(l)}_x/M_2) \end{array} \;\; \Longleftrightarrow \;\; \begin{array}{c} \text{Fourier space}\ (N^{(k)}_x,\, N^{(k)}_z/M_2,\, N^{(k)}_y/M_1) \end{array}$$2, where $$\begin{array}{c} \text{Physical space}\ (N^{(l)}_y,\, N^{(l)}_z/M_1,\, N^{(l)}_x/M_2) \end{array} \;\; \Longleftrightarrow \;\; \begin{array}{c} \text{Fourier space}\ (N^{(k)}_x,\, N^{(k)}_z/M_2,\, N^{(k)}_y/M_1) \end{array}$$3 is the reachable subspace of $$\begin{array}{c} \text{Physical space}\ (N^{(l)}_y,\, N^{(l)}_z/M_1,\, N^{(l)}_x/M_2) \end{array} \;\; \Longleftrightarrow \;\; \begin{array}{c} \text{Fourier space}\ (N^{(k)}_x,\, N^{(k)}_z/M_2,\, N^{(k)}_y/M_1) \end{array}$$4 (Ferrante et al., 2012). The decomposition then yields explicit parameterizations of optimal state and control trajectories through a finite-dimensional linear endpoint system.

These examples illustrate a recurring distinction. In PDE codes and MRI, pencil decomposition is primarily a computational device for scaling or signal separation. In operator theory, tensor decomposition, and optimal control, it is a structural transformation whose purpose is to isolate spectrum, constraints, reachable components, or canonical blocks. A plausible implication is that the phrase is best interpreted not by its literal wording alone, but by the ambient object being decomposed: a distributed grid, a sampled time series, or a parameter-dependent linear system.