Biquadratic Optimization with Orthogonality Constraints

• The paper introduces a canonical model for biquadratic optimization with orthogonality constraints, demonstrating that every non-optimal critical point is a strict saddle.
• It explains how reformulations on the Stiefel manifold and Riemannian geometry enable a clear landscape structure and efficient parameterizations.
• Algorithmic frameworks like quasi-Newton updates and constraint-dissolving methods reduce computational costs while ensuring convergence in joint diagonalization and PCA applications.

Biquadratic optimization with orthogonality constraints studies nonconvex matrix problems in which the decision variable is constrained by an orthogonality relation—most commonly $X^\top X=I$—while the objective is biquadratic, quartic, or becomes quartic after reformulation. A canonical model is the homogeneous quadratic program on the Stiefel manifold,

$$\max_{\bm Q \in \mathrm{St}(d,K)} \operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right),$$

which recent landscape theory treats explicitly as a matrix biquadratic form over orthonormal columns (Wang et al., 25 Jun 2026). Closely related formulations arise in approximate joint diagonalization, where orthogonality-constrained reformulation and low-rank compression yield the JADOC objective

$$\min_{B\in O(N)} \frac{1}{2K}\sum_{k=1}^K\sum_{i=1}^N \log\!\left(\lambda+\sum_{j=1}^S (BL_k)_{ij}^2\right),$$

and in the classical Stiefel-constrained problem

$\min_{X\in {\rm St}(m,n)} {\rm tr}(X^TAXB)$

that underlies much of the modern convergence theory (Vlaming et al., 2021, Liu et al., 2015).

1. Canonical problem classes

The most explicit recent formulation is the quadratic program with orthogonality constraints (QPOC),

$$\max_{\bm Q \in \mathrm{St}(d,K)} F(\bm Q):=\operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right), \qquad \mathrm{St}(d,K)=\{\bm Q\in\mathbb R^{d\times K}:\bm Q^T\bm Q=\bm I_K\}.$$

If $\bm q_1,\dots,\bm q_K$ are the columns of $\bm Q$ and $\bm B=(b_{ij})$, then

$$\operatorname{Tr}(\bm Q^T\bm A\bm Q\bm B)=\sum_{i,j=1}^K b_{ij}\,\bm q_i^T\bm A\,\bm q_j.$$

The literature characterizes this as a canonical biquadratic optimization problem because the objective is quartic in the entries of $\bm Q$, bilinear in pairs of columns, and quadratic in $$\max_{\bm Q \in \mathrm{St}(d,K)} \operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right),$$0; when $$\max_{\bm Q \in \mathrm{St}(d,K)} \operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right),$$1 and $$\max_{\bm Q \in \mathrm{St}(d,K)} \operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right),$$2 are diagonal, it becomes a weighted matrix biquadratic form over orthonormal columns. This model includes eigenvalue problems, generalized eigenvalue problems, and orthogonal relaxations of quadratic assignment problems (Wang et al., 25 Jun 2026).

A second major class is approximate joint diagonalization of symmetric positive semidefinite matrices $$\max_{\bm Q \in \mathrm{St}(d,K)} \operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right),$$3. The starting criterion recalled in JADOC is

$$\max_{\bm Q \in \mathrm{St}(d,K)} \operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right),$$4

which measures the gap between the determinant of the diagonal part and the full determinant. Under orthogonality, the problem is reformulated into an orthogonality-preserving optimization over a reduced representation, and the paper emphasizes that the resulting objective is a biquadratic optimization problem with orthogonality constraints in the sense that the transformed matrices depend quadratically on $$\max_{\bm Q \in \mathrm{St}(d,K)} \operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right),$$5, while the criterion depends on logarithms of sums of squares of their entries (Vlaming et al., 2021).

The broader orthogonality-constrained optimization class is

$$\max_{\bm Q \in \mathrm{St}(d,K)} \operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right),$$6

with smooth, possibly nonconvex $$\max_{\bm Q \in \mathrm{St}(d,K)} \operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right),$$7. Recent geometric work makes explicit that this framework applies to quadratic and biquadratic forms as special cases, even when the objective itself is not singled out as biquadratic in the paper title or abstract (Goyens et al., 21 Jul 2025).

2. Orthogonality as manifold geometry

The standard feasible set is the Stiefel manifold,

$$\max_{\bm Q \in \mathrm{St}(d,K)} \operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right),$$8

with tangent space

$$\max_{\bm Q \in \mathrm{St}(d,K)} \operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right),$$9

Equivalent representations include $$\min_{B\in O(N)} \frac{1}{2K}\sum_{k=1}^K\sum_{i=1}^N \log\!\left(\lambda+\sum_{j=1}^S (BL_k)_{ij}^2\right),$$0, which is central for update parameterizations based on skew-symmetric generators (Goyens et al., 21 Jul 2025).

In JADOC, the decisive reformulation is to impose

$$\min_{B\in O(N)} \frac{1}{2K}\sum_{k=1}^K\sum_{i=1}^N \log\!\left(\lambda+\sum_{j=1}^S (BL_k)_{ij}^2\right),$$1

Then $$\min_{B\in O(N)} \frac{1}{2K}\sum_{k=1}^K\sum_{i=1}^N \log\!\left(\lambda+\sum_{j=1}^S (BL_k)_{ij}^2\right),$$2, so the second log-determinant term becomes independent of $$\min_{B\in O(N)} \frac{1}{2K}\sum_{k=1}^K\sum_{i=1}^N \log\!\left(\lambda+\sum_{j=1}^S (BL_k)_{ij}^2\right),$$3 and drops from the objective. This restricts the search to the orthogonal group, keeps iterates orthonormal throughout, and enables the update

$$\min_{B\in O(N)} \frac{1}{2K}\sum_{k=1}^K\sum_{i=1}^N \log\!\left(\lambda+\sum_{j=1}^S (BL_k)_{ij}^2\right),$$4

where $$\min_{B\in O(N)} \frac{1}{2K}\sum_{k=1}^K\sum_{i=1}^N \log\!\left(\lambda+\sum_{j=1}^S (BL_k)_{ij}^2\right),$$5 is strictly lower triangular and $$\min_{B\in O(N)} \frac{1}{2K}\sum_{k=1}^K\sum_{i=1}^N \log\!\left(\lambda+\sum_{j=1}^S (BL_k)_{ij}^2\right),$$6 is skew-symmetric. The matrix exponential therefore guarantees orthogonality by construction (Vlaming et al., 2021).

Recent work also generalizes orthogonality to constraints of the form

$$\min_{B\in O(N)} \frac{1}{2K}\sum_{k=1}^K\sum_{i=1}^N \log\!\left(\lambda+\sum_{j=1}^S (BL_k)_{ij}^2\right),$$7

with $$\min_{B\in O(N)} \frac{1}{2K}\sum_{k=1}^K\sum_{i=1}^N \log\!\left(\lambda+\sum_{j=1}^S (BL_k)_{ij}^2\right),$$8 and $$\min_{B\in O(N)} \frac{1}{2K}\sum_{k=1}^K\sum_{i=1}^N \log\!\left(\lambda+\sum_{j=1}^S (BL_k)_{ij}^2\right),$$9 a linear self-adjoint mapping. The feasible set

$\min_{X\in {\rm St}(m,n)} {\rm tr}(X^TAXB)$0

is shown to be a closed embedded submanifold. Its tangent and normal spaces are

$\min_{X\in {\rm St}(m,n)} {\rm tr}(X^TAXB)$1

and this framework covers the Stiefel, generalized Stiefel, symplectic Stiefel, indefinite Stiefel, hyperbolic, and third-order tensor Stiefel manifolds (Zhang et al., 29 Apr 2026).

A complementary geometric development is the introduction of a family of metrics

$\min_{X\in {\rm St}(m,n)} {\rm tr}(X^TAXB)$2

which extends the classical $\min_{X\in {\rm St}(m,n)} {\rm tr}(X^TAXB)$3-metric from the Stiefel manifold to full-rank matrices. On the Stiefel manifold, $\min_{X\in {\rm St}(m,n)} {\rm tr}(X^TAXB)$4 recovers the Euclidean metric and $\min_{X\in {\rm St}(m,n)} {\rm tr}(X^TAXB)$5 the canonical metric. This metric family underlies a landing framework in which a tangent term decreases the objective and a normal term reduces infeasibility (Goyens et al., 21 Jul 2025).

3. Critical points and landscape structure

For QPOC, diagonal reduction is the basic structural simplification. If

$\min_{X\in {\rm St}(m,n)} {\rm tr}(X^TAXB)$6

then $\min_{X\in {\rm St}(m,n)} {\rm tr}(X^TAXB)$7 preserves the Stiefel constraint and transforms the objective into diagonal coordinates. One may therefore assume without loss of generality that $\min_{X\in {\rm St}(m,n)} {\rm tr}(X^TAXB)$8 and $\min_{X\in {\rm St}(m,n)} {\rm tr}(X^TAXB)$9 are diagonal with eigenvalues sorted in nonincreasing order (Wang et al., 25 Jun 2026).

The Riemannian critical set of QPOC is characterized by

$$\max_{\bm Q \in \mathrm{St}(d,K)} F(\bm Q):=\operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right), \qquad \mathrm{St}(d,K)=\{\bm Q\in\mathbb R^{d\times K}:\bm Q^T\bm Q=\bm I_K\}.$$0

When $$\max_{\bm Q \in \mathrm{St}(d,K)} F(\bm Q):=\operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right), \qquad \mathrm{St}(d,K)=\{\bm Q\in\mathbb R^{d\times K}:\bm Q^T\bm Q=\bm I_K\}.$$1 is nonsingular, every critical point admits a closed-form description involving orthogonal blocks, a basis-selection matrix, and a permutation. When $$\max_{\bm Q \in \mathrm{St}(d,K)} F(\bm Q):=\operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right), \qquad \mathrm{St}(d,K)=\{\bm Q\in\mathbb R^{d\times K}:\bm Q^T\bm Q=\bm I_K\}.$$2 is singular but nonzero, the critical set decomposes into a reduced QPOC on the nonzero part of $$\max_{\bm Q \in \mathrm{St}(d,K)} F(\bm Q):=\operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right), \qquad \mathrm{St}(d,K)=\{\bm Q\in\mathbb R^{d\times K}:\bm Q^T\bm Q=\bm I_K\}.$$3 together with free orthogonal degrees of freedom in the null part. This explicit parameterization is the basis for the paper’s full landscape classification (Wang et al., 25 Jun 2026).

The central theorem is that every critical point of QPOC is either a global maximizer or a strict saddle point. Global optimality is characterized by an exact matching condition $$\max_{\bm Q \in \mathrm{St}(d,K)} F(\bm Q):=\operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right), \qquad \mathrm{St}(d,K)=\{\bm Q\in\mathbb R^{d\times K}:\bm Q^T\bm Q=\bm I_K\}.$$4, where $$\max_{\bm Q \in \mathrm{St}(d,K)} F(\bm Q):=\operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right), \qquad \mathrm{St}(d,K)=\{\bm Q\in\mathbb R^{d\times K}:\bm Q^T\bm Q=\bm I_K\}.$$5 is the blockwise sorted selection of diagonal entries of $$\max_{\bm Q \in \mathrm{St}(d,K)} F(\bm Q):=\operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right), \qquad \mathrm{St}(d,K)=\{\bm Q\in\mathbb R^{d\times K}:\bm Q^T\bm Q=\bm I_K\}.$$6 at the critical point and $$\max_{\bm Q \in \mathrm{St}(d,K)} F(\bm Q):=\operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right), \qquad \mathrm{St}(d,K)=\{\bm Q\in\mathbb R^{d\times K}:\bm Q^T\bm Q=\bm I_K\}.$$7 is the ideal assignment pairing positive eigenvalues of $$\max_{\bm Q \in \mathrm{St}(d,K)} F(\bm Q):=\operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right), \qquad \mathrm{St}(d,K)=\{\bm Q\in\mathbb R^{d\times K}:\bm Q^T\bm Q=\bm I_K\}.$$8 with the largest diagonal entries of $$\max_{\bm Q \in \mathrm{St}(d,K)} F(\bm Q):=\operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right), \qquad \mathrm{St}(d,K)=\{\bm Q\in\mathbb R^{d\times K}:\bm Q^T\bm Q=\bm I_K\}.$$9, and negative eigenvalues with the smallest. If this condition fails, the paper constructs a tangent direction of positive curvature using Givens rotations, thereby proving strict saddleness of every non-optimal critical point (Wang et al., 25 Jun 2026).

The same paper identifies the population version of heteroscedastic probabilistic PCA (HePPCA) as a special instance of QPOC. In that model, the global maximizers are exactly

$\bm q_1,\dots,\bm q_K$0

so the true principal subspace is recovered up to column-wise sign flips. Moreover, the population objective is locally geodesically strongly concave near every global maximizer, and the sample problem inherits the same benign landscape with high probability when the sample size is sufficiently large (Wang et al., 25 Jun 2026).

4. Reformulations and algorithmic frameworks

The JADOC algorithm combines orthogonality, low-rank approximation, and quasi-Newton updates. Each $\bm q_1,\dots,\bm q_K$1 is approximated by a rank-$\bm q_1,\dots,\bm q_K$2 factorization $\bm q_1,\dots,\bm q_K$3, regularized as $\bm q_1,\dots,\bm q_K$4 with

$\bm q_1,\dots,\bm q_K$5

and the recommended rank is

$\bm q_1,\dots,\bm q_K$6

The search direction is obtained from a diagonal quasi-Newton approximation,

$\bm q_1,\dots,\bm q_K$7

with Hessian entries clipped below $\bm q_1,\dots,\bm q_K$8 to avoid numerical instability. The line search is accelerated by linearizing the update

$\bm q_1,\dots,\bm q_K$9

so each candidate $\bm Q$0 costs only $\bm Q$1 (Vlaming et al., 2021).

Landing methods provide a different strategy. The update

$\bm Q$2

combines a tangent descent direction for the objective with a normal correction term for infeasibility. For the $\bm Q$3-metric family, the projection onto the tangent space and the constrained gradient admit explicit formulas, which makes the method applicable to any smooth orthogonality-constrained objective, including biquadratic ones (Goyens et al., 21 Jul 2025).

Constraint-dissolving methods remove orthogonality constraints by embedding them into an unconstrained penalty. In the generalized setting, the constraint-dissolving mapping is

$\bm Q$4

and the unconstrained objective is

$\bm Q$5

This construction satisfies $\bm Q$6 on the feasible set and yields equivalence with the original constrained problem in the sense of first-order and second-order stationary points, locally for sufficiently large $\bm Q$7 (Zhang et al., 29 Apr 2026).

Decentralized orthogonality-constrained optimization introduces communication and sampling costs in addition to nonconvex geometry. VRSGT avoids explicit Riemannian projections and retractions by using an augmented Lagrangian estimate together with the practical local direction

$\bm Q$8

where $\bm Q$9 coincides with the Riemannian gradient on the Stiefel manifold. Gradient tracking variables then combine consensus and variance reduction over the network (Wang et al., 2022).

When nonnegativity is added to orthogonality, the feasible set acquires a strong combinatorial structure: every feasible matrix has at most one nonzero entry per row. The Support-Set algorithm exploits this by fixing a support pattern and solving the proximal linearization subproblem in closed form. The resulting iterates remain strictly feasible throughout and the support pattern is updated by relocating nonzero entries only when the current support is insufficiently descent-producing (Wang et al., 5 Nov 2025).

5. Complexity and convergence properties

For the classical Stiefel-constrained quadratic problem, the foundational rate result is an explicit Łojasiewicz inequality with exponent $\bm B=(b_{ij})$0: $\bm B=(b_{ij})$1 near a critical point $\bm B=(b_{ij})$2. The proof relies on a local error bound,

$\bm B=(b_{ij})$3

and on an explicit description of the critical set after simultaneous diagonalization of $\bm B=(b_{ij})$4 and $\bm B=(b_{ij})$5. Combined with standard descent conditions, this yields linear convergence of a broad class of retracted Armijo line-search methods, including methods based on the polar decomposition retraction, QR retraction, Cayley transform, and the Riemannian exponential map (Liu et al., 2015).

JADOC’s principal efficiency claim is a reduction from the $\bm B=(b_{ij})$6 per-iteration cost typical of existing joint diagonalization methods to $\bm B=(b_{ij})$7 per quasi-Newton iteration. The $\bm B=(b_{ij})$8 work is moved into one-time low-rank preprocessing, after which the shared orthogonal update and the compressed factors allow joint updates of all $\bm B=(b_{ij})$9 terms at $$\operatorname{Tr}(\bm Q^T\bm A\bm Q\bm B)=\sum_{i,j=1}^K b_{ij}\,\bm q_i^T\bm A\,\bm q_j.$$0 cost. Empirically, JADOC is dramatically faster than qndiag, which is itself much faster than JADE, while JADE performs best in diagonalization quality, JADOC is very close behind, and qndiag is substantially worse (Vlaming et al., 2021).

The generalized penalty approach also emphasizes per-iteration savings. Computing $$\operatorname{Tr}(\bm Q^T\bm A\bm Q\bm B)=\sum_{i,j=1}^K b_{ij}\,\bm q_i^T\bm A\,\bm q_j.$$1 for GOCDF costs

$$\operatorname{Tr}(\bm Q^T\bm A\bm Q\bm B)=\sum_{i,j=1}^K b_{ij}\,\bm q_i^T\bm A\,\bm q_j.$$2

whereas computing the Riemannian gradient for a general $$\operatorname{Tr}(\bm Q^T\bm A\bm Q\bm B)=\sum_{i,j=1}^K b_{ij}\,\bm q_i^T\bm A\,\bm q_j.$$3 may require solving a least-squares problem with cost as large as

$$\operatorname{Tr}(\bm Q^T\bm A\bm Q\bm B)=\sum_{i,j=1}^K b_{ij}\,\bm q_i^T\bm A\,\bm q_j.$$4

On the indefinite Stiefel manifold, the paper reports $$\operatorname{Tr}(\bm Q^T\bm A\bm Q\bm B)=\sum_{i,j=1}^K b_{ij}\,\bm q_i^T\bm A\,\bm q_j.$$5 for the Riemannian gradient and $$\operatorname{Tr}(\bm Q^T\bm A\bm Q\bm B)=\sum_{i,j=1}^K b_{ij}\,\bm q_i^T\bm A\,\bm q_j.$$6 for the Euclidean gradient via CDF (Zhang et al., 29 Apr 2026).

In decentralized optimization, VRSGT attains

$$\operatorname{Tr}(\bm Q^T\bm A\bm Q\bm B)=\sum_{i,j=1}^K b_{ij}\,\bm q_i^T\bm A\,\bm q_j.$$7

when $$\operatorname{Tr}(\bm Q^T\bm A\bm Q\bm B)=\sum_{i,j=1}^K b_{ij}\,\bm q_i^T\bm A\,\bm q_j.$$8, with communication rounds $$\operatorname{Tr}(\bm Q^T\bm A\bm Q\bm B)=\sum_{i,j=1}^K b_{ij}\,\bm q_i^T\bm A\,\bm q_j.$$9 and total sampled gradients $\bm Q$0 to reach an $\bm Q$1-stationary point. The paper contrasts this with DRGTA, which has sample complexity $\bm Q$2, and DRSGD, which has sample complexity $\bm Q$3 (Wang et al., 2022).

For orthogonality plus nonnegativity, Support-Set has subproblems of cost $\bm Q$4, overall per-iteration cost at most $\bm Q$5, global convergence to a first-order stationary point, and iteration complexity $\bm Q$6 for reaching an $\bm Q$7-approximate first-order stationary point (Wang et al., 5 Nov 2025).

The benign-landscape theorem for QPOC adds a complementary perspective: because every nonoptimal critical point is a strict saddle and prior work provides a Łojasiewicz exponent $\bm Q$8, retraction-based methods are expected to avoid saddles from random initialization, converge to global maximizers, and do so at a local linear rate (Wang et al., 25 Jun 2026).

6. Applications, variants, and scope

Approximate joint diagonalization remains a central application domain. JAD is used in estimation of common principal components, estimation of multiple variance components, and blind signal separation. The JADOC reformulation shows that these tasks can be attacked by orthogonality-constrained quasi-Newton optimization with low-rank compression, rather than by repeated full-matrix sweep operations (Vlaming et al., 2021).

HePPCA provides a statistically structured instance of QPOC. The model

$\bm Q$9

leads to a nonconvex maximum-likelihood problem on the Stiefel manifold, and the population version becomes a diagonal QPOC in coordinates aligned with $$\max_{\bm Q \in \mathrm{St}(d,K)} \operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right),$$00. This makes landscape results directly relevant to latent-subspace recovery under heteroscedastic noise (Wang et al., 25 Jun 2026).

Structured variants broaden the scope of the field. For decentralized PCA and decentralized DPCP, orthogonality constraints interact with consensus constraints and stochastic sampling; VRSGT reports promising performance in a real-world autonomous driving application based on KITTI point clouds (Wang et al., 2022). For nonnegative PCA, clustering, community detection, orthogonal nonnegative matrix factorization, spectral clustering, nonnegative Laplacian embedding, copositivity checking, and quadratic assignment problems, the combination of nonnegativity and orthogonality induces row-wise one-sparsity that can be exploited algorithmically (Wang et al., 5 Nov 2025).

The recent literature also distinguishes between explicit biquadratic models and broader orthogonality-constrained smooth optimization frameworks. The landing framework does not formulate a specific biquadratic objective, but it applies directly to any smooth $$\max_{\bm Q \in \mathrm{St}(d,K)} \operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right),$$01, including quadratic or biquadratic forms, under $$\max_{\bm Q \in \mathrm{St}(d,K)} \operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right),$$02 (Goyens et al., 21 Jul 2025). The generalized constraint-dissolving approach likewise targets a wider class of quadratic matrix orthogonality constraints of the form $$\max_{\bm Q \in \mathrm{St}(d,K)} \operatorname{Tr}\!\left(\bm Q^T \bm A \bm Q \bm B\right),$$03 (Zhang et al., 29 Apr 2026). A plausible implication is that biquadratic optimization with orthogonality constraints is now studied both as a specific algebraic model and as a special case of a larger manifold-optimization and penalty-method toolkit.