The Normal Procrustes Problem: A Riemannian Optimization Approach

Abstract: For given $m \times n$ data matrices $X, Y$, we investigate the Normal Procrustes Problem---the least squares optimization problem that aims to minimize $|AX-Y|_F<sup>2$, where $A$ is constrained to be a normal $m \times m$ matrix. As far as the author of this article is aware, no other method that attempts to solve the Normal Procrustes Problem exists in the literature; we thus propose what is, to our knowledge, the first such method. We, furthermore, adapt our approach to address the Real Normal Procrustes Problem, where $A$ must be real. In our treatment of these problems, we first reduce our complex and real objective functions to be purely optimizable over the Riemannian manifolds of the unitary and real orthogonal matrices, respectively. This reduction enables us to apply techniques in Riemannian manifold optimization to approximate solutions to both. The Closest Normal Matrix and Real Closest Normal Matrix Problems are both special cases of their respective Procrustes Problems and have been previously studied in the literature. Our approach thus recovers a novel Riemannian optimization method for approximating solutions to both these problems. We further numerically test the performance of our method across all such problems (including against previously developed algorithms on the Closest Normal Matrix Problems) and obtain competitive residuals and favorable scaling in wall-clock time.

• The paper introduces the first general computational method for the Normal Procrustes Problem, reducing it to optimization over unitary or orthogonal matrices with explicit Riemannian gradients and Hessians.
• The complex solver recovered planted normal solutions with residuals below 10⁻¹⁰ in every trial, matched established methods on the Closest Normal Matrix Problem, and scaled favorably for larger matrices.
• The real formulation is substantially more difficult because block pairings create nonsmoothness and suboptimal critical points, with only 45.8% of recovery trials reaching residuals below 10⁻¹⁰.

Problem statement and motivation

For given $X, Y \in M_{m\times n}(\mathbb{C})$ with $m \ge n$, the Normal Procrustes Problem (NPP) asks for a normal matrix $A \in \mathcal{N} = \{B : B^*B = BB^*\}$ minimizing $\|AX - Y\|_F^2$. The paper addresses this problem because no prior method exists in the literature specifically targeting it; the author claims this is the first such method. The motivation stems from the $\mathcal{P}$-Targeting Problem—determining when $AX = Y$ exactly—which remains unsolved for normal $A$, largely because the underlying Normal Matrix Completion Problem lacks a closed-form solution. Rather than pursuing exact conditions, the paper approximates $A^\star$ via Riemannian manifold optimization.

The paper also treats the Real NPP ($A \in \mathcal{N}_R$, real data), and notes that setting $X = I_m$ recovers the well-studied Closest Normal Matrix Problem (CNP), enabling comparison against Ruhe's Jacobi-type algorithm and Guglielmi–Scalone's two-level gradient system.

Reduction to unitary optimization

Exploiting the spectral decomposition $m \ge n$0 with $m \ge n$1 and unitary invariance of the Frobenius norm, the objective decouples across eigenvalues:

$m \ge n$2

For fixed $m \ge n$3, each $m \ge n$4 (or zero when $m \ge n$5), where $m \ge n$6 and $m \ge n$7. The main reduction theorem shows the NPP is equivalent to maximizing

$m \ge n$8

over $m \ge n$9, recovering $A \in \mathcal{N} = \{B : B^*B = BB^*\}$0. Geometrically, $A \in \mathcal{N} = \{B : B^*B = BB^*\}$1 is a weighted sum of squared cosines of Hermitian angles between corresponding rows of $A \in \mathcal{N} = \{B : B^*B = BB^*\}$2 and $A \in \mathcal{N} = \{B : B^*B = BB^*\}$3. Crucially, the paper demonstrates via an explicit geodesic example that $A \in \mathcal{N} = \{B : B^*B = BB^*\}$4 is not geodesically concave, so Riemannian gradient methods cannot guarantee global convergence—a limitation that shapes all subsequent numerical claims.

Riemannian gradient and Hessian: complex case

On the open dense set where all $A \in \mathcal{N} = \{B : B^*B = BB^*\}$5 (all of $A \in \mathcal{N} = \{B : B^*B = BB^*\}$6 when $A \in \mathcal{N} = \{B : B^*B = BB^*\}$7), the Euclidean gradient is

$A \in \mathcal{N} = \{B : B^*B = BB^*\}$8

projected onto $A \in \mathcal{N} = \{B : B^*B = BB^*\}$9 via the skew-Hermitian part operator. The Riemannian Hessian follows Boumal's standard formula after differentiating the Euclidean gradient, requiring derivatives $\|AX - Y\|_F^2$0 obtained by the quotient rule. Manopt's checkgradient/checkhessian confirm slopes of approximately 2 and 3 respectively across four rank regimes at $\|AX - Y\|_F^2$1, with tangency, linearity, and symmetry residuals vanishing to machine precision.

Numerical results: complex case

Using Manopt's trustregions solver with random initialization on $\|AX - Y\|_F^2$2, experiments over $\|AX - Y\|_F^2$3 and $\|AX - Y\|_F^2$4 show several notable findings:

• Known-optimum recovery: in every trial where a normal $\|AX - Y\|_F^2$5 with $\|AX - Y\|_F^2$6 was planted, the solver attained final residual below $\|AX - Y\|_F^2$7; exact recovery of $\|AX - Y\|_F^2$8 itself occurred whenever $\|AX - Y\|_F^2$9, while smaller $\mathcal{P}$0 yielded equally sufficient but distinct solutions due to freedom on $\mathcal{P}$1.
• Non-uniqueness: for small $\mathcal{P}$2, different starts achieve identical residuals but distinct matrices (spread up to ~127% of $\mathcal{P}$3).
• Suboptimal stationary points: occasional cases exhibit differing residuals across starts, directly evidencing the failure of geodesic concavity; tightening the gradient tolerance to $\mathcal{P}$4 did not collapse these differences.
• Anomalous scaling: solve time peaks near $\mathcal{P}$5 rather than at $\mathcal{P}$6; the paper concedes the reason for this asymmetry "remains unclear."

Comparison on the CNP

Setting $\mathcal{P}$7 reduces the objective to $\mathcal{P}$8. Against Ruhe's algorithm (reimplemented from pseudocode, since his code is unavailable) and Guglielmi–Scalone's method, all three methods agree on final residuals for benchmark matrices from both papers—with one exception: on the Frank matrix $\mathcal{P}$9, the paper's method and Guglielmi–Scalone's attain marginally lower residuals than the reimplemented Ruhe's. Notably, the paper corrects a claim in Guglielmi–Scalone that Ruhe's algorithm "is proved to determine the global optimum": Ruhe only proves convergence to a stationary point with local optimality conditions.

A consistent structural advantage emerges in the normality measure $AX = Y$0: Guglielmi–Scalone's iterates exceed $AX = Y$1 in every trial (normality only in the limit), whereas the proposed method stays below $AX = Y$2. On scaling tests up to $AX = Y$3, Ruhe's method is fastest for $AX = Y$4, but the crossover occurs at $AX = Y$5, beyond which the trust-region solver is fastest in every instance (e.g., 352 s vs. 1387 s vs. 5290 s at $AX = Y$6). A log–log fit gives an iteration-count exponent of roughly 0.7 in $AX = Y$7, so observed wall-clock cost scales between $AX = Y$8 and $AX = Y$9.

Real case: reduction and additional pathologies

Real normal matrices admit a real quasidiagonal form $A$0 with $A$1 and rotation-form $A$2 blocks. For fixed $A$3, each block is optimized independently, choosing between a rotation block $A$4 and a diagonal block by comparing residuals. This yields a reduced objective $A$5 involving a pointwise maximum over the two branch types, maximized over $A$6.

Two sources of non-smoothness arise: the measure-zero loci where denominators vanish (handled as in the complex case) and, uniquely, a tie locus $A$7 where the two branches of the maximum coincide. Random initialization avoids $A$8 with probability one, though the paper does not prove that solver trajectories avoid it.

Gradient and Hessian checks again pass. However, the known-optimum recovery experiment reveals substantially worse behavior than the complex case: only 45.8% of trials achieved residual below $A$9 (62.5% below $A^\star$0). The paper attributes this to increased nonconcavity from loss of permutation invariance: because $A^\star$1 has block structure, coordinate pairings matter. An explicit $A^\star$2 example constructs a suboptimal critical point $A^\star$3 where the Riemannian gradient vanishes and the Hessian is negative semidefinite (eigenvalues $A^\star$4 and $A^\star$5), yet the residual is 20 versus the global optimum of 0. With $A^\star$6 possible pairings for even $A^\star$7, many such traps may exist, though the paper does not quantify how many are actually suboptimal critical points.

Comparison on the Real CNP

Against Guglielmi–Scalone on real benchmarks, results diverge more than in the complex case: their method achieves lower residuals on three instances, the proposed method wins on $A^\star$8, and residuals differ by 3–20% elsewhere. Overall, Guglielmi–Scalone obtained lower residuals on 6 of 9 tested instances, so the paper explicitly declines to claim superiority in solution quality; its advantages are speed and exact satisfaction of normality. The paper suggests that multiple random starts, affordable given the speed advantage, could close the residual gap—an empirical question left open.

Limitations and open questions

Several limitations are stated plainly in the paper. Global optimality is never guaranteed: $A^\star$9 and $A \in \mathcal{N}_R$0 are not geodesically concave, whether the infimum of the NPP is always attained is not determined, and smoothness fails on measure-zero sets whose avoidance along solver trajectories is unproven. The anomalous solve-time peak near $A \in \mathcal{N}_R$1 is unexplained. In the real case, suboptimal termination is frequent, and the efficacy of remedies (multiple starts, geometry-aware initialization avoiding degeneracies for rank-deficient $A \in \mathcal{N}_R$2) is deferred to future work. Extension to the Normal Matrix Completion Problem with prescribed entry patterns, which Guglielmi–Scalone's framework accommodates, is also left open. Finally, comparisons against Ruhe's algorithm rely on the author's own reimplementation, which may not match a tuned original.

Conclusion

The paper establishes the first computational method for the Normal Procrustes Problem in full generality, reducing it to smooth (almost everywhere) optimization over $A \in \mathcal{N}_R$3 or $A \in \mathcal{N}_R$4 with explicit Riemannian gradients and Hessians suitable for second-order trust-region methods. Empirically, the complex-case solver is highly reliable—recovering planted global optima to residual below $A \in \mathcal{N}_R$5 in every trial—and matches established CNP algorithms' solution quality while scaling favorably beyond $A \in \mathcal{N}_R$6. The real case is markedly harder, with frequent suboptimal termination traced to permutation-asymmetric block structure inducing spurious critical points; improving initialization and characterizing these traps remain open problems.