Conceptio › Archive › arXiv CS

arXiv CSopen access

Deflation-Free Optimal Scoring

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

neural-networks

machine learning, deep learning, neural networks

Deflation-Free Sparse Optimal Scoring Sharmin Afroz1 and Brendan Ames2

arXiv:2604.25664v1 [stat.ML] 28 Apr 2026

Department of Computational & Informational Sciences, Stillman College, Tuscaloosa, AL, USA, [email protected] 2 School of Mathematical Sciences, University of Southampton, Southampton, UK, [email protected] April 29, 2026 Abstract Sparse Optimal Scoring (SOS) reformulates linear discriminant analysis to enable feature selection through elastic net regularization, making it well-suited for high-dimensional settings where the number of features exceeds observations. Most existing SOS methods use deflation-based strategies that compute discriminant vectors sequentially, which can propagate errors and produce suboptimal solutions. We propose a novel approach that estimates all discriminant vectors simultaneously under an explicit global orthogonality constraint, which we call Deflation-Free Sparse Optimal Scoring (DFSOS). DFSOS combines Bregman iteration with orthogonality-constrained optimization, decomposing the problem into tractable subproblems for scoring vectors, discriminant vectors, and orthogonality enforcement. We establish convergence to stationary points of the augmented Lagrangian under mild conditions. Extensive experiments using synthetic data and real-world time series data demonstrate that DFSOS achieves classification accuracy comparable to or better than existing deflation-based methods. These results indicate that deflation-free approaches offer a robust and effective framework for sparse discriminant analysis in high-dimensional problems.

Introduction

Linear Discriminant Analysis (LDA) is a classical supervised learning technique widely used for classification, dimension reduction, and data visualization. Given a dataset with n observations, p features, and K classes, LDA seeks a low-dimensional linear subspace in which class separation is maximized. Under the multivariate Gaussian model, LDA assumes that observations from class i follow a normal distribution N (µi , Σ) with a common covariance matrix across classes. Classification is performed by maximizing posterior class probabilities, which leads to linear decision boundaries. This Bayesian formulation yields discriminant functions of the form 1 δi (x) = xT Σ−1 µi − µTi Σ−1 µi + log πi . 2 where πi denotes the prior probability of class i. We will focus on an equivalent framework for LDA based on optimal scoring [16, 15]. Here, LDA is reformulated as a linear regression problem with numerical scores assigned to class labels. Under appropriate constraints, this optimal scoring approach is known to coincide with the Bayesian formulation (see [17, Section 4.3] and [26]). Fisher proposed an alternative formula of LDA based on maximizing the Rayleigh quotient of the between-class and within-class scatter matrices of the training data. This reduces LDA to a generalized eigenvalue problem. Unfortunately, when the number of features p is larger than the number of observations n, the sample covariance matrix and sample within-class scatter matrices are singular; in this case, the change of basis/variables needed to perform LDA is impossible to compute. This shortcoming of LDA inspired penalized or sparse discriminant analysis (SDA) approaches which introduce regularization, particularly sparsity-inducing penalties, into the discriminant analysis framework so that only a subset of informative features is used for classification [15, 38, 30, 12, 27, 14, 7, 36, 9]. Many existing sparse discriminant analysis methods find multiple discriminant directions using a deflationary process. After each discriminant direction is computed, the next one is forced to be orthogonal to the previous ones. While this step-by-step process is relatively easy to implement, it has important 1

limitations. Error in the early directions can carry over and affect the later ones, which can hurt both interpretability and classification performance. In addition, this sequential procedure makes parallel computation difficult and hides the overall structure of the optimization problem. These limitations motivate the development of deflation-free methods. Instead of computing discriminant vectors one at a time, these methods estimate all directions at once while directly enforcing orthogonality between them. This method maintains sparsity while avoiding the numerical and statistical problems caused by sequential deflation.

1.1

Contribution

We propose a novel formulation for sparse discriminant analysis that jointly estimates multiple discriminant vectors under a global orthogonality constraint, eliminating the need for sequential deflation and avoiding error propagation inherent in deflation-based approaches. This approach is based on a modification of the sparse optimal scoring approach considered in [9, 2]. Specifically, we reformulate the deflation-based sparse optimal scoring problem using a global orthogonality constraint. We then propose and analyse a block-coordinate descent method for solution of this orthogonality-constrained optimization problem. We establish conditions for convergence of this heuristic and observe via empirical trials that it is competitive in terms of classification performance and computational complexity.

1.2

Outline

In Section 2, we review sparse optimal scoring and existing deflation-based block coordinate descent approaches, including their optimization properties and limitations. In Section 3, we introduce the proposed deflation-free sparse optimal scoring formulation and present the DFSOS algorithm, including detailed derivations of each subproblem, implementation details, and convergence analysis. In Section 4, we evaluate the empirical performance of DFSOS on synthetic Gaussian data and real time-series classification benchmarks, comparing accuracy, sparsity, and computational efficiency with existing methods. Finally, Section 5 concludes the paper with a discussion of implications, limitations, and directions for future research.

Sparse Optimal Scoring

2.1

The Optimal Scoring Problem

As mentioned previously, the optimal scoring problem seeks to recast the classification problem as a regression problem by using a sequence of scorings to turn the categorical class-labels into quantitative variables. Here, we formally discuss the specialization of this framework to linear discriminant analysis. Let X ∈ Rn×p denote a data matrix with rows encoding n observations of p-dimensional data vectors and let Y be an n × K indicator matrix for the K classes. It is defined as Yij = 1 if the ith observation belongs to the jth class, Yij = 0 otherwise. Let θ i ∈ RK denote the ith scoring vector and β i ∈ Rp be the ith discriminant vector. The optimal scoring problem seeks to jointly optimize the scoring-discriminant vector (θ i , β i ) so that the sum of squared differences between class labels encoded by θ i and the linear regression model encoded by β i is minimized: min

∥Y θ i − Xβ i ∥2

s.t.

T 1 T n θi Y Y θi = 1 T T θi Y Y θℓ = 0

β i ,θ i

(1) ∀ℓ < i.

We assume that the scoring vector θ i is normalized and orthogonal to the previously computed i − 1 scoring vectors to make the problem well-posed. This process is inherently deflationary as we must compute scoring-discriminant vector pairs one by one. LDA has been observed to yield accurate classifiers when applied to data with a small number of predictor variables, particularly for Gaussian data when each class has a shared population covariance matrix. However, it can overfit the data in problems with large numbers of highly correlated variables. Conversely, when the class boundaries are nonlinear it can underfit the data [15]. Moreover, when the number of features p is much larger than the number of observations n i.e., p > n, the sample covariance matrix Σ, proportional to X T X, does not have full rank. To address these issues, Clemmensen et al. [9]

introduce regularization in the form of an elastic net penalty function [38]. This yields the sparse or penalized optimal scoring problem (SOS): min

∥Y θ i − Xβ i ∥2 + γβ Ti Ωβ i + λ∥β i ∥1

s.t.

T 1 T n θi Y Y θi = 1 T T θi Y Y θℓ = 0

β i ,θ i

(2) ∀ℓ < i.

It utilizes both the lasso penalty (ℓ1 -norm), to induce sparsity, and the ridge penalty (ℓ2 -norm), to encourage selection of clusters of correlated variables [38]. More details on the use of the naive elastic net penalty can also be found in [17, Section 18.4]. Note that the generalized elastic net penalty differs from the naive elastic net only in the use of the ℓ2 -norm with respect to a different basis.

2.2

A Block Coordinate Descent Approach for Sparse Optimal Scoring

Clemmensen et al. [9] propose a heuristic for solving (1) using block coordinate descent with respect to θ i and β i . To update θ i , we fix β i and solve the problem minθi s.t.

∥Y θ i − Xβ i ∥2 T 1 T n θi Y Y θi = 1 T T θ i Y Y θ ℓ = 0 ∀ℓ < i.

(3)

Although this is a non-convex quadratic program in θ i , (3) admits a closed-form solution, given in the following theorem. Theorem 1. The problem (3) has optimal solution θ new = s(I − Qi QTi D)D −1 Y T Xβ i i

(4)

where D = n1 Y T Y , Qi is the K × i matrix with columns consisting of the previous i − 1 scoring vectors θ 1 , θ 2 , ..., θ i−1 and the all-ones vector 1 ∈ RK , and s is a proportionality constant ensuring that θ Ti Dθ i = 1. In particular, θ i is given by w = (I − Qi QTi D)D −1 Y T Xβ i ,

θi = √

w wT Dw

(5)

After updating θ i , we fix θ i and solve the β-subproblem = argmin ∥Y θ i − Xβ∥2 + γβ T Ωβ + λ∥β∥1 . β new i

(6)

This is a convex optimization problem in β and can be solved efficiently using iterative methods. Clemmensen et al. propose the use of least angle regression (LARS) [11, 38, 28] to solve (6). Alternately, Atkins et al. [2] consider solution of (6) using the accelerated proximal gradient method and the alternating direction of method of multipliers. In either case, we solve (6) and then use this new value of β to update θ and so on. We repeat this process until the iterates (θ i , β i ) converge. After finding the optimal scoring-discriminant vector pairs, classification is performed using nearest centroid classification on the projection of the data onto the space spanned by the discriminant vectors, given by [Xβ 1 . . . Xβ q ] where q < K [9]. This iterative approach for updating θ and β is summarized in Algorithm 1. The following theorem establishes that the sequence of function values for SOS given by Algorithm 1 is convergent; see [2, Theorem 2.4]. Theorem 2. Suppose that the sequence of iterates {(θ t , β t )}∞ t=0 is generated by Algorithm 1. Then the defined by F (θ, β) := ∥Y θ−Xβ∥2 +γβ T Ωβ+λ∥β∥1 sequence of objective function values {F (θ t , β t )}∞ t=0 is convergent. It can also be established that every convergent subsequence of {(θ t , β t )}∞ t=1 converges to a stationary point of (2), see [2, Theorem 2.5]. tℓ tℓ ∞ Theorem 3. Let {(θ t , β t )}∞ t=1 be the sequence of points generated by Algorithm 1. Suppose that {(θ , β )}ℓ=1 t t ∞ ∗ ∗ ∗ ∗ is a convergent subsequence of {(θ , β )}t=1 with limit (θ , β ). Then (θ , β ) is a stationary point of (2): (θ ∗ , β ∗ ) is feasible for (2) and there exists ψ ∗ ∈ R and ν ∗ ∈ Rj−1 such that 0 ∈ ∂L(θ ∗ , β ∗ , ψ ∗ , ν ∗ ), where ∂L(θ, β, ψ, ν) denotes the subdifferential of the Lagrangian function L with respect to the primal variables (θ, β).

Algorithm 1 Sparse Discriminant Analysis [9, 2] 1: Given stopping tolerance ϵ and maximum number of iterations N . ∗ ∗ ∗ ∗ ∗ ∗ 2: return scoring-discriminant vector pairs (θ 1 , β 1 ), (θ 2 , β 2 ), . . . , (θ K−1 , β K−1 ). 3: for j = 1, 2, . . . , K − 1 do 4: Initialize θ 0j :

θ 0j = 5: 6: 7: 8:

1 I − Qj QTj Y T Y (Y T Y )−1 z, n

√

θ 0j =

nθ 0j

∥Y θ 0j ∥

Calculate the jth scoring and discriminant vector pair (θ ∗j , β ∗j ): for i = 0, 1, 2 . . . N do Update β ij as the solution of (6) with θ = θ ij . Update θ i+1 by j w=

I−

1 Q QT Y T Y n j j

√

(Y T Y )−1 Y T Xβ ij ,

θ i+1 =

nw ; ∥Y w∥

Converged if the residual between consecutive iterates is smaller than tolerance:

10:

if max

∥θ i+1 −θ ij ∥ ∥β i+1 −β ij ∥ j j ∥θ i+1 ∥ j

∥β i+1 ∥ j

< ϵ then

11: The algorithm has converged. 12: Increment j and break 13: end if 14: end for 15: end for

Deflation-Free Sparse Optimal Scoring

Most existing SOS algorithms rely on deflation-based strategies. Deflation-based methods impose orthogonality implicitly and locally, rather than enforcing a global orthogonality constraint across all discriminant vectors. As a result, the final set of discriminant directions may be sensitive to initialization, ordering, and numerical errors, and may fail to optimally capture class separation in a joint sense. These challenges become more pronounced in high-dimensional settings ( p ≫ n), where sparse regularization already introduces nonconvexity and instability. To overcome these limitations, deflation-free methods reformulate SOS as a single orthogonality-constrained optimization problem, allowing all discriminant and scoring vectors to be computed simultaneously. By enforcing orthogonality at the matrix level and avoiding sequential deflation, deflation-free approaches eliminate error accumulation, improve numerical stability, and yield more coherent discriminant subspaces. This motivates the development of algorithms such as Deflation-Free Sparse Optimal Scoring (DFSOS), which leverage splitting techniques and Bregman iteration to efficiently handle orthogonality constraints while preserving sparsity and scalability in high-dimensional settings. We consider the following variant of the sparse optimal scoring problem, which we call deflation-free sparse optimal scoring: min

J(θ, β) := ∥Y θ − Xβ∥2F + γβ T β + λ∥β∥1

s.t.

T 1 T n θ Y Y θ = I.

θ∈RK×q ,β∈Rp×q

(7)

Note that (7) requires minimization of a quadratic function in (θ, β) over the manifold of K ×q orthogonal matrices (with respect to the inner product ⟨θ1 , θ2 ⟩ = n1 θ1 T Y T Y θ2 ). Note further that we substitute Ω = I in the ridge regression term in (2); our problem and algorithm can be modified to allow generalized Tikhonov regression, but we focus on this specialization for the sake of model and algorithmic simplicity. This problem is non-convex. Moreover, optimization over the manifold of orthogonal matrices is NP-hard; see [21, 22] and the set of reductions given in [19, Section 1].

3.1

A Splitting Method for Deflation-Free Sparse Optimal Scoring

Due to this worst-case intractability, much recent research has focused on the design of heuristics for solving orthogonally-constrained optimization problems [1, 20, 35, 13, 23, 25, 18, 8, 4, 3, 37]. In particular, we will leverage a method recently proposed by Lai and Osher [20] to approximately solve (7). This approach uses matrix splitting and Bregman iteration, to generate a sequence of approximate solutions of (7) from the solution of easier to solve subproblems. We introduce the primal variable P = √1n Y θ. Then the DFSOS problem can be written as arg min

J(θ, β)

θ∈RK×q ,β∈Rp×q

P = √1n Y θ,

s.t.

(8)

P T P = I.

Using Bregman iteration the optimal solution of (8) can be found iteratively using the update formulas ρ (θ k , β k , P k ) = arg min J(θ k , β k ) + ∥LX − P k−1 + B k−1 ∥2F , 2 θ, β, P ∈Rn

s.t.

P T P = I,

B k = B k−1 + LX k − P k ,

(9) (10)

where L := √1n Y . We will approximately solve the subproblem (9) using block-coordinate descent, where we minimize with respect to each primal variable θ, β, and P individually with the others fixed.

3.2

The θ sub-problem

We first fix β, P , B and describe the process for updating θ. In this case, θ new is the solution of min θ

s.t.

q X

∥Y θ i − Xβ i ∥2 + γ∥β∥2 + λ∥β∥1 +

i=1 T T

θ Y Y 1 = 0,

ρ 1 √ Y θi − P i + B i 2 n

i = 1, 2, ...., q.

, the ith column of θ or ith This problem is separable with respect to the columns of θ. Thus, θ new i scoring vector, is the minimizer of min θ

s.t.

q X

∥Y θ i − Xβ i ∥2 + γ∥β i ∥2 + λ∥β i ∥1 +

i=1 θ Ti Y T Y 1 = 0.

ρ 1 √ Y θi − P i + B i 2 n

The objective function of Equation (11) is a convex quadratic in θ i : ρ 1 T T 2 T T T T T T √ θ Y Y θi − θ Y (P i − B i ) + C F (θ i ) = θ i Y Y θ i − 2θ i Y Xβ i + 2 n i n i

(11)

(12)

where C represents all the constant terms depending only on β, P and B. Therefore, the sub-problem for θ i can be reformulated as: T T ρ min 1 + 2n θ i Y Y θ i − 2θ Ti Y T Xβ i + 2√ρ n (P i − B i ) + C θi (13) s.t. θ Ti Y T Y 1 = 0. This is a strongly convex quadratic program in θ i . Therefore, the unique solution of the Karush-KuhnTucker conditions is the unique minimizer. Enforcing the KKT conditions yields the following linear system for (θ i , vi ), where vi is the Lagrange multiplier corresponding to the equality constraint:      T 2Y (Xβ i + 2√ρ n (P i − B i )) 2I + nρ Y T Y −Y T Y 1 θi     =  (14) .  T T v 1 Y Y 0 0 i Thus, we update θ i as the solution of linear system (14) using the current values of β i , P i , B i for each i = 1, 2, ...., q. Note that the coefficient matrix in (14) depends only on the parameter ρ and the matrix Y T Y ; we can pre-compute a Cholesky factorization and use it to update each θ i until ρ is changed. This update scheme is also amenable to parallelization as we can update all θ i simultaneously. 5

Algorithm 2 DFSOS Algorithm 0

1: Initialize: Given regularization parameter λ > 0 and initial θ , β , set B

= 0, P 0 = LX, where

L = √1n Y . 2: while “not converged” do 3: (θ k , β k ) = arg min J(θ k−1 , β k−1 ) + ρ2 ∥ √1n Y θ k−1 − P k−1 + B k−1 ∥2F .

Compute the SVD Lθ k + B k−1 = U DV T . 5: P k = U I n×m V T . 6: B k = B k−1 + Lθ k − P k 7: k ←k+1 8: end while 4:

3.3

Updating β

We update β new as the solution of the following unconstrained optimization problem min β

q X

∥Y θ i − Xβ i ∥2 + γ∥β i ∥2 + λ∥β i ∥1 ,

(15)

i=1

with θi fixed, where βi denotes the ith column of β. This objective function is separable in βi . Each βi can be updated independently using the iterative methods outlined in Section 2.2. In the experimental analysis found later in this manuscript, we will focus on the use of the accelerated proximal gradient method proposed by Atkins et al. [2]; however, any efficient iterative method for non-smooth convex minimization may be used. Again, we can parallelize this update and compute all new βi simultaneously.

3.4

Updating P

Finally, we update P new by minimizing the augmented Lagrangian of (8) with respect to P with θ and β fixed. The following theorem gives a closed-form solution for P new . Theorem 4 ([20, Theorem 2.1]). The constrained quadratic problem: 1 ∥Q − Y ∥2F , Q∈Rn×m 2

Q∗ = arg min

s.t.

QT Q = I

(16)

has an analytical solution Q∗ = U I n×m V T , where U ∈ Rn×n , V ∈ Rm×m are orthogonal matrices and D ∈ Rn×m is a diagonal matrix satisfying the singular value decomposition (SVD) Y = U DV T . Moreover, if rank(Y ) = m, then Q∗ = Y Ṽ D̃

−1/2

Ṽ ,

where V ∈ Rm×m is an orthogonal matrix and D̃ ∈ Rm×m is a diagonal matrix satisfying the SVD T Y T Y = Ṽ D̃ Ṽ . Note that P new = argmin ∥P − (Lθ + B)∥2F

s.t.

PTP = I

P ∈Rk×q

for fixed iterates (θ, β, B). Theorem 4 suggests an algorithm for updating P : compute the singular value decomposition Lθ + B = U DV T ; then let P new = U I k×q V T .

3.5

(17)

A Heuristic for Deflation-Free Optimal Scoring

Each iteration of our proposed algorithm DFSOS updates (θ, β, P ) using (14), (15), and (17), followed by the approximate dual ascent step (10). We summarize the steps of the DFSOS Algorithm in Alg. 2.

3.6

Convergence of DFSOS

Convergence of our algorithm follows directly from Corollary 2 given in [34, Section 5]. We have the following theorem: Theorem 5. Let {(θ k , β k )}∞ k=1 be the sequence of points generated by Algorithm 2. Then DFSOS (Algorithm 2) is convergent for sufficiently large choice of r. In this case, each set of iterates {(θ k , β k )}∞ k=1 is bounded, has at least one limit point, and each limit point (θ ∗ , β ∗ ) is a stationary point of the augmented Lagrangian L(θ, β, P ), i.e., 0 ∈ ∂L(θ ∗ , β ∗ , P ∗ ). Proof. The objective function J(θ, β) is strongly convex and differentiable. For l1 -constrained optimization, Bregman iteration coincides with the general ADMM Algorithm. It is easy to see that both ADMM and Bregman Iteration are bounded and the augmented Lagrangian is lower bounded at a critical point [34, 24]. Moreover, each set of the sequence {(θ k , β k )}∞ k=1 generated by Algorithm 2 is a bounded set and J(θ, β) is lower bounded on the feasible set. This shows that Algorithm 2 and the problem (8) satisfy the hypothesis of Corollary 2 [34]. This establishes the stated convergence of Algorithm 2.

3.7

Practical Implementation Details

We conclude this section with a discussion of practical details guiding usage of our heuristic for deflationfree sparse optimal scoring. 3.7.1

Range of Regularization Parameters

We want to choose λ so we that are guaranteed a non-trivial solution in the β subproblem. If we set λ = 0, the β-subproblem (15) becomes min β

q X

∥Y θ i − Xβ i ∥2 + γ∥β i ∥2 ,

(18)

i=1

where θ is fixed. Expanding gives the equivalent problem min β

q X

β Ti X T Xβ i − 2β Ti X T Y θ i + γ∥β i ∥2 .

(19)

i=1

Applying stationarity and taking partial derivatives with respect to β shows that β ∗ is a solution of the linear system 2 X T X + γI β ∗ = 2X T Y θ (20) This linear system has unique solution −1 β ∗ = X T X + γI X T Y θi .

(21)

Substituting β ∗ into the original objective function (19) gives the minimum value q X

∗ ∗T T T T T β ∗T X X + γI β − 2β X Y θ + θ Y Y θ + λ∥β ∥ , i i 1 i i i i i

(22)

i=1

where we have θ Ti Y T Y θ i = n by feasibility of θ. Now, rearranging we have a nontrivial β solution if β ∗ has value at most nq. Therefore, we need q X

T β ∗T X T X + γI β ∗i − 2β ∗T i i X Y θ i + λ∥β i ∥1 ≤ 0.

(23)

i=1

Thus, it is sufficient to choose Pq λ≤

T T T T i=1 β i (X X + γI)β i − 2β i X Y θ i

∥β∥1 7

(24)

Algorithm 3 Pseudocode for initializing θ 1: Input random θ ∈ Rk×q . 2: for i = 1, 2, . . . do 3: D = n1 (Y T Y ) 4: t = u − Qj (QTj Du) 5: a = tT Dt √ 6: θ i = t/ a 7: end for to ensure there exists at least one solution β ∗ with objective value strictly less than the value of the trivial solution β = 0. A more flexible approach is to assign a different regularization term λi for each column of β: min β

q X

∥Y θ i − Xβ i ∥2 + γ∥β i ∥2 + λi ∥β i ∥1 .

(25)

i=1

In this case, we can ensure a non-trivial solution by choosing λ1 , λ2 , . . . , λq to satisfy λi ≤

β Ti (X T X + γI)β i − 2β Ti X T Y θ i ∥β i ∥1

for all i = 1, 2, . . . , q. 3.7.2

Initialization

We randomly initialize θ and apply an orthogonalization process based on modification of the GramSchmidt process to ensure that our initial iterate θ satisfies θ T Y T Y θ = nI. Details of this process are given in Algorithm 3. On the other hand, we calculate initial β from the initial orthogonal θ and the data matrix X using the Sherman-Morrison-Woodbury process. Lemma 6 (Sherman-Morrison-Woodbury). If A = E + U V , then A−1 = (E + U V )−1 = E −1 − E −1 U (I + V E −1 U )−1 V E −1 .

(26)

After temporarily setting the ℓ1 -penalty parameter λ = 0, (21) is equivalent to min β

q X

β Ti X T Xβ i − 2β Ti X T Y θ ∗i + γ∥β i ∥2 .

(27)

i=1

Here θ ∗i is the orthogonal initial θ. Applying stationarity and taking partial derivatives of (27) with respect to β i gives optimum β ∗ such that 2 X T X + γI β ∗i = 2X T Y θ ∗i . (28) Applying the Sherman-Morrison-Woodbury formula (26), shows that −1 β ∗i = X T X + γI X T Y θ ∗i −1 ! 12 T 1 1 T I− X I + XX X X T Y θ ∗i = γ γ γ ! −1 1 1 T 1 T = M− X I + XX (XM ) , γ γ γ where M = X T Y θ ∗i . This process is summarized in Alg. 4.

Algorithm 4 Pseudocode for β ∗ 1: Input X, θ , Y , γ ∗ T 2: Calculate B = X Y θ i 3: Compute XB T 4: Compute I + γ1 XX 5: Using Cholesky factorization solve

I + γ1 XX T V = XB

1. RT R = I + γ1 XX T 2. Solve RT RV = XB for RV ⇒ temp = RV = RT \(XB) 3. Solve V = R\temp T 6: Compute β = γ1 B − γ1 X V . Algorithm 5 Update of penalty parameter ρ 0

1: Input η < 1, σ > 1, initial ρ, initial orthogonal primal variable P 0 = Lθ . 2: Initialize v0 = 2∥P 0 ∥2F = 2q.

3: At iteration k: Compute v = ∥P − Lθ∥2F 4: if v < ηvk then 5: vk+1 = v 6: else 7: ρk+1 = σρk 8: vk+1 = vk 9: end if

3.7.3

An Alternative Factorization

The matrix Y T Y is a diagonal matrix with ith diagonal entry equal to the number of observations ni belonging to class i. We can perform the matrix-splitting given in (8) using the matrix √  n1 √  n2 1    (29) L= √   ..  n . √ nK . This factorization can be more computationally efficient than that outlined above as the matrix L is k × k, rather than n × n. Moreover, products involving L only require scaling according to the diagonal entries of L. 3.7.4

The Augmented Lagrangian Parameter

We dynamically update the augmented Lagrangian parameter r each iteration. If the residual ∥P k − Lθ k ∥2F following iteration k does not decrease a sufficient amount, then we increase r to accelerate convergence of the Bregman iteration scheme as suggested in [6, Section 3.1]. Details for this update process are given in Alg. 5.

Empirical Analysis

We next compare the empirical performance of our deflation-free optimal scoring methods with existing methods for SOS, as well as standard classification methods, for a variety of data.

4.1

Gaussian Data

We first investigate the performance of our algorithm for classifying Gaussian data. We performed the following experiment 10 times for each (K, r)-pair for K ∈ {3, 6} and r ∈ {0.1, 0.5, 0.9}. In each experiment, we sample p-dimensional vectors from K multivariate normal distributions for p = 1000. For each dataset, we apply sparse optimal scoring to perform dimension reduction, and validate this 9

transformation by performing nearest centroid classification following projection onto the span of the learned discriminant vectors. 4.1.1

Distribution of Gaussian Data

We generate training data corresponding to the ith class, i = 1, 2, ..., K, by sampling ntrain = 100 observations from the multivariate normal distribution with mean µi ∈ Rp satisfying ( 0.7 if 100(i − 1) < j ≤ 100i [µi ]j = 0 otherwise for all j = 1, 2, ..., K and covariance matrix Σ ∈ Rp×p with all features correlated defined by ( r, if i ̸= j, Σij = 1, if i = j for fixed scalar r ∈ [0, 1). In each experiment, we sampled ntest = 1000 testing observations from each class. For each combination of (K, r), we produced 10 training-validation pairs of data. 4.1.2

Methods for Comparison and Choice of Parameters

We apply the DFSOS method proposed in Section 3 to fit discriminant vectors using each training data 1/2 set. We consider two variants based on solving (8) with L = √1n Y (DFSOS-1 ) and L = √1n Y T Y (DFSOS-2 ). We use the accelerated proximal gradient heuristic proposed in [2] to solve (10) and update β in each version of DFSOS. We compare our deflation-free approach with the accelerated proximal gradient (APG) and alternating direction method of multipliers (ADMM) heuristics for sparse optimal scoring proposed in [2]. We use 5-fold cross validation to choose ℓ1 -penalty parameter λ in DFSOS, APG, and ADMM from 7 possible choices on the exponential grid −3 −2 −1 2 , 2 , 2 , 1, 2, 22 , 23 × λmax , where λmax is chosen as in Section 3.7.1. We also compare our classification results using support vector machine (SVM) and k-nearest neighbour (KNN) classifiers for k = 1, 5, 10. DFSOS requires choice of stopping criteria for the inner optimization of β and the overall optimization iteration. We terminate the inner optimization if a 10−5 suboptimal solution is found or a maximum of 50 iterations are performed; we stop the β-update proximal gradient step after a 10−4 suboptimal solution is found or after 100 iterations. We terminate the outer loop after a maximum of 500 iterations or a 10−4 suboptimal solution is found. We use scaling parameters η = 0.25, σ = 2, and initial ρ = 5 to dynamically update the augmented Lagrangian parameter in DFSOS as described in Section 3.7.4. We choose initial β and θ as in Section 3.7.2. The ADMM and APG heuristics require choice of the regularization parameters Ω, γ, and λ. In each experiment, we set Ω = I and γ = 10−1 . We choose λ by 5-fold cross validation from an exponential grid equivalent with that used for DFSOS. These heuristics require stopping criteria for inner optimization of β and the outer optimization loop. We terminate the inner loop after a maximum of 100 iterations or if a 10−4 suboptimal solution is found. We terminate the outer loop if a 10−4 suboptimal solution is found or a maximum of 500 iterations have been performed. The augmented Lagrangian parameter for ADMM was set at µ = 2 for all experiments. 4.1.3

Computational Environments

All experiments were performed with Matlab R2023a [32] using a single serial node of the IRIDIS 5 High Performance Computing Facility at the University of Southampton. Matlab implementations of the DFSOS and ASDA methods are shared at Dr. Ames’s Github repository1 . We used default settings for the support vector machine and nearest neighbours functions from Matlab’s Statistics and Machine Learning Toolbox [33, 31].

4.2

Discussion of Experiments involving Gaussian Data

Summary statistics for our experiments involving Gaussian data can be found in Table 1. 1 https://github.com/bpames/Deflation-Free-SOS

Data

Measure

DFSOS-1

DFSOS-2

APG

ADMM

SVM

1-NN

5-NN

10-NN

K=3 r = 0.1

Accuracy Run-time (s) Cardinality

1.0 (0.0) 4.325 (1.107) 0.119 (0.0)

0.999 (0.0) 3.948 (0.512) 0.113 (0.0)

0.999 (0.0) 2.462 (0.127) 0.118 (0.0)

1.0 (0.0) 2.626 (0.25) 0.166 (0.023)

1.0 (0.0) 0.174 (0.014) –

0.92 (0.0) 0.021 (0.0) –

0.978 (0.0) 0.029 (0.001) –

0.988 (0.0) 0.022 (0.0) –

K=3 r = 0.5

Accuracy Run-time (s) Cardinality

1.0 (0.0) 5.084 (1.373) 0.084 (0.0)

1.0 (0.0) 4.783 (0.635) 0.085 (0.0)

0.999 (0.0) 2.87 (0.068) 0.097 (0.0)

0.998 (0.0) 3.099 (0.265) 0.093 (0.0)

1.0 (0.0) 0.17 (0.012) –

0.98 (0.0) 0.021 (0.0) –

0.981 (0.0) 0.027 (0.001) –

0.973 (0.0) 0.021 (0.0) –

K=3 r = 0.9

Accuracy Run-time (s) Cardinality

1.0 (0.0) 5.022 (1.783) 0.097 (0.001)

1.0 (0.0) 4.695 (0.656) 0.113 (0.001)

0.588 (0.002) 2.446 (0.059) 0.012 (0.0)

0.614 (0.007) 2.554 (0.234) 0.012 (0.0)

1.0 (0.0) 0.167 (0.012) –

0.991 (0.0) 0.022 (0.0) –

0.976 (0.0) 0.028 (0.001) –

0.962 (0.0) 0.021 (0.0) –

K=6 r = 0.1

Accuracy Run-time (s) Cardinality

0.996 (0.0) 12.095 (0.561) 0.168 (0.0)

0.994 (0.0) 11.283 (0.482) 0.17 (0.0)

0.992 (0.0) 12.761 (0.177) 0.141 (0.0)

0.993 (0.0) 12.737 (0.545) 0.144 (0.0)

1.0 (0.0) 0.66 (0.015) –

0.855 (0.0) 0.036 (0.0) –

0.949 (0.0) 0.045 (0.001) –

0.978 (0.0) 0.034 (0.0) –

K=6 r = 0.5

Accuracy Run-time (s) Cardinality

0.999 (0.0) 15.44 (0.876) 0.191 (0.001)

0.999 (0.0) 14.614 (0.536) 0.203 (0.001)

0.942 (0.001) 16.253 (0.276) 0.069 (0.0)

0.951 (0.001) 16.179 (1.285) 0.069 (0.0)

1.0 (0.0) 0.655 (0.014) –

0.969 (0.0) 0.038 (0.0) –

0.977 (0.0) 0.046 (0.001) –

0.972 (0.0) 0.035 (0.0) –

K=6 r = 0.9

Accuracy Run-time (s) Cardinality

1.0 (0.0) 17.329 (0.992) 0.123 (0.002)

1.0 (0.0) 16.439 (0.668) 0.124 (0.002)

0.942 (0.003) 14.416 (0.401) 0.094 (0.0)

0.943 (0.003) 14.24 (0.935) 0.095 (0.0)

1.0 (0.0) 0.64 (0.014) –

0.99 (0.0) 0.036 (0.0) –

0.975 (0.0) 0.044 (0.001) –

0.966 (0.0) 0.032 (0.0) –

Table 1: Out-of-sample prediction accuracy, run-time in seconds, and fraction of non-zero features for discriminant vectors and classifiers for a selection of Gaussian data constructed as described in Section 4.1.1. Data is reported in the format mean (variance) over 10 experiments for each (K, r)-pair. 4.2.1

Comparison of Classifiers

Our proposed deflation-free methods consistently provide more accurate classifiers than the deflationary approaches based on the ASDA algorithm, with comparable run-times. Run-times for all four methods were within a few seconds of each other in all experiments, while discriminant vectors computed using DFSOS tended to be slightly more dense than those given by ASDA. Moreover, both DFSOS methods matched or exceeded the accuracy of the SVM and KNN classifiers in all experiments. We must note that SVM and KNN classifiers required much less time to fit, as they do not require any training of regularization parameters. This phenomenon was starkest when the data was highly correlated, e.g., when r = 0.9. For example, when r = 0.9 and K = 3, the ASDA classifiers both struggled to identify meaningful discriminant directions. In each case, the ASDA method converged to a pair of discriminant vectors with a modest number of relatively small non-zero entries. To illustrate, we consider a random set of Gaussian data constructed as described in Section 4.1 with K = 3 and r = 0.9. The discriminant vectors β given by each method for these data are visualised in Fig. 1. We can see that the discriminant vectors found using ASDA have much smaller magnitude compared to those found using DFSOS. Moreover, the support of the pairs of discriminant vectors found by ASDA differ from those found by DFSOS, which have many more nonzero entries; we should also note that while the discriminant vectors found by DFSOS-1 and DFSOS-2 appear to differ, this is due to rotational symmetry in the paired scoring vectors and the dimension reduction and classifiers for the two DFSOS approaches behave very similarly. This shrinkage of ASDA discriminant vectors to 0 leads to significant degradation of classification performance. We visualise decision boundaries of the nearest centroid classifier following projection onto the span of discriminant vectors in Fig. 3. We can clearly observe that there is minimal separation between classes in the space spanned by the discriminant vectors found by ASDA (see Fig. 3c and Fig. 3d). This is largely avoided when we use DFSOS. Here, observations in each class are well-separated in the discriminant vector space. We compared this behaviour with that for a set of Gaussian data constructed as described in Section 4.1 with K = 3 and r = 0.5. Visualisations of discriminant vectors and classification boundaries are given in Fig. 2 and Fig. 4 respectively. We can see that all found methods give roughly the same discriminant vectors, decision boundaries, and predictions, up to rotational symmetry. However, even in this case, we see that the test data is better separated in the span of the discriminant vectors given by DFSOS. 4.2.2

Hypothesis Tests for Gaussian Data

For each (K, r)-pair, we performed one-sided Wilcoxon tests to compare pairs of classification methods in terms of prediction error, run-times, and cardinality. For measure f , we test the null hypothesis H0 : f (i) = f (j) against the one-sided alternative Ha : f (i) < f (j) for each pair of methods i and j. The

'LVFULPLQDQW9HFWRUV*DXVVLDQV(k = 3, r = 0.9)')626

(b) DFSOS-2 'LVFULPLQDQW9HFWRUV*DXVVLDQV(k = 3, r = 0.9)$'00

(a) DFSOS-1 'LVFULPLQDQW9HFWRUV*DXVVLDQV(k = 3, r = 0.9)$3* '9

(d) ASDA-ADMM

'9 '9

Figure 1: Visualisation of discriminant vectors returned by each of DFSOS-1, DFSOS-2, and ASDA with APG and ADMM subproblem solvers. Note that the nonzero entries of the discriminant vectors found by ASDA have much small magnitude than that of the discriminant vectors found by each variant of DFSOS; this is indicated by very pale colours in Fig. 1c and Fig. 1d.

(b) DFSOS-2

(a) DFSOS-1

(d) ASDA-ADMM

Figure 2: Visualisation of discriminant vectors returned by each of DFSOS-1, DFSOS-2, and ASDA with APG and ADMM subproblem solvers. Note that the nonzero entries of the discriminant vectors found by ASDA, unlike in Fig. 1, have similar magnitude to the discriminant vectors found by each variant of DFSOS.

'HFLVLRQ%RXQGDULHVDQG3UHGLFWLRQV*DXVVLDQV(k = 3, r = 0.9)')626

&ODVV

'LVFULPLQDQW9HFWRU

(a) DFSOS-1

(b) DFSOS-2

'HFLVLRQ%RXQGDULHVDQG3UHGLFWLRQV*DXVVLDQV(k = 3, r = 0.9)$3*

'HFLVLRQ%RXQGDULHVDQG3UHGLFWLRQV*DXVVLDQV(k = 3, r = 0.9)$'00

&ODVV

'LVFULPLQDQW9HFWRU

&ODVV

'LVFULPLQDQW9HFWRU

(d) ASDA-ADMM

Figure 3: Nearest centroid classification decision boundaries and predictions for testing data in the span of each discriminant vector pair for K = 3, r = 0.9. Actual class label is indicated by colour and predicted class label is indicated by shape. Note that test data is not separable when using the discriminant vectors given by ASDA.

'LVFULPLQDQW9HFWRU

&ODVV

'LVFULPLQDQW9HFWRU

&ODVV

(a) DFSOS-1

(b) DFSOS-2 &ODVV

'LVFULPLQDQW9HFWRU

&ODVV

'LVFULPLQDQW9HFWRU

(d) ASDA-ADMM

Figure 4: Nearest centroid classification decision boundaries and predictions for testing data in the span of each discriminant vector pair for K = 3, r = 0.5. Actual class label is indicated by colour and predicted class label is indicated by shape. Note that test data is more separable when using the discriminant vectors given by DFSOS.

690

$3 *

69 0

') 62 6

(a) Classification Accuracy

(b) Run-times

$'00

$3 *

') 62 6

')626

$3*

$'00

$3*

') 62 6

$3 *

') 62 6

1XPEHURIH0UHMHFWLRQV:LOFR[RQ&DUGLQDOLW\

')626

$3*

$'00

') 62 6

1XPEHURIH0UHMHFWLRQV:LOFR[RQ5XQ7LPH

')626 ')626

1XPEHURIH0UHMHFWLRQV:LOFR[RQ$FFXUDF\

Figure 5: Number of (K, r)-pairs for which we observe statistically significant improvement in classification accuracy, computational efficiency, and cardinality of discriminants when using row method compared to column method. That is, the (i, j)-entry is the number of (K, r)-pair where we reject the null hypothesis H0 : there is no difference between method i and method j using the Wilcoxon test with one-sided alternative hypothesis Ha : method i is better than method j. number of times we reject H0 for each (i, j)-pair using an α = 0.05 level of significance for each measure is given in Fig. 5. We observe the following: • Classifiers calculated using the two variants of DFSOS were significantly more accurate than all methods except for SVM in most trials. DFSOS classifiers were more accurate than all KNN classifiers for all (K, r); DFSOS-1 was significantly more accurate than both ASDA methods for all but one (K, r)-pair, and DFSOS-2 was significantly more accurate than both ASDA methods for all but two (K, r)-pairs. We should note that we observe similar phenomena using other accuracy metrics such as recall, precision, and f -score. • DFSOS tended to find denser discriminant vectors than ASDA; ASDA had significantly sparse discriminant vectors for four (K, r) pairs. • DFSOS required, on average, more time to converge than the ASDA methods, albeit by a small amount of actual time (1-2 seconds in most trials). ASDA had significantly less run-time for four (K, r)-pairs. On the other hand, DFSOS-2 required less computation on average than DFSOS-1 for all (K, r)-pairs. From these hypothesis tests, we can conclude that DFSOS provides significantly better classification accuracy than ASDA, at a modest cost of increased run-time and cardinality of solutions. Moreover, this classification performance matches the state of the art for linearly separable data given by support vector machines, with increased interpretability due to the sparsity of discriminant vectors. 4.2.3

Consistency of Predictions

Although we have observed an improvement in classification performance using DFSOS (and SVM) over ASDA in these experiments, we should caution that the actual classification performance of all of the observed methods is remarkably similar. To illustrate this phenomenon, we calculate the cosine similarity for the out-of-sample predictions made by each classifier. The cosine similarity between the vector of predictions for each pair of classifiers, averaged over the 10 problem instances sampled for each (K, r)pair, is given in Fig. 6. The diagonal entries of Fig. 6 indicate the average cosine similarity between all pairs of prediction vectors made by each method across all 10 Gaussian observations. The cosine similarities for the class membership predictions summarized in Fig. 6 are almost all nearly equal to 1, even when the Wilcoxon test suggests a statistically significant difference in classification accuracy. This suggests that when there is a difference in classification prediction, it is due to a disagreement in a small number of assigned class labels, with all other predictions agreeing for the two classifiers. This can be observed in Fig. 4. In this setting, we have (K, r) = (3, 0.5) and our hypothesis tests suggest that DFSOS classifiers provided a statistically significant improvement in classification performance over ASDA classifiers. However, the decision boundaries for the classifiers are similar and 15

690

$' 00

0 69

1 1

') 6

(d) (K, r) = (6, 0.1)

(e) (K, r) = (6, 0.5)

1 1

$'00

1 1

$' 00

$3*

') 6

')626 ')626

$3*

$'00

690

$YHUDJH&RVLQH6LPLODULW\6\QWK

') 6

6 ')

1 1

0 69

$YHUDJH&RVLQH6LPLODULW\6\QWK

6 ')

62 ')

(b) (K, r) = (3, 0.5)

')626 ')626

$YHUDJH&RVLQH6LPLODULW\6\QWK

(a) (K, r) = (3, 0.1)

') 6

$' 00

') 6

')626 ')626

$3*

$'00

690

')626 ')626

$3*

$'00

690

') 6

$'00

$3*

')626 ')626

$3*

$'00

690

$YHUDJH&RVLQH6LPLODULW\6\QWK

')626 ')626

$YHUDJH&RVLQH6LPLODULW\6\QWK

(f) (K, r) = (6, 0.9)

Figure 6: Average cosine similarity between out-of-sample predictions for each pair of methods. The diagonal entries indicate the average cosine similarity between all pairs of prediction vectors made by each method across all 10 Gaussian observations. the predicted class labels agree for all four DFSOS and ASDA classifiers except for a relatively small number of misclassified observations near the decision boundary for the two ASDA classifiers.

4.3

Empirical Analysis of Time-Series Data

We also evaluated performance of DFSOS using a subset of data from the UC-Riverside classification data repository [10]. We limited our analysis to data containing more predictor variables than observations (p > n) and at least K = 3 classes; this yielded a collection of 15 data sets. We used DFSOS, APG, and ADMM to obtain q = K − 1 sparse discriminant vectors and performed nearest-centroid classification after projection onto the subspace spanned by these discriminant vectors. We compare classification accuracy of these methods with each other and standard classifiers SVM, KNN as before. For each data set we train each classifier using 10 random initializations and compare average run-time and classification accuracy across these 10 trials. 4.3.1

Experimental Parameters

As before, we use 5-fold cross validation to choose the ℓ1 -penalty λ from the exponential grid {2−4 , 2−3 , 2−2 , 2−1 , 1} × λmax in DFSOS, APG, and APG. We terminate the iterative update of β in DFSOS, APG, and ADMM after 50 iterations or a 10−5 suboptimal solution is found and we terminate the outer loop after 50 iterations or a 10−4 suboptimal solution is found. We use η = 0.25, σ = 2, and initial ρ = 5 to update the augmented Lagrangian parameter ρ in DFSOS. We use augmented Lagrangian parameter µ = 2 in ADMM.

Data

Measure

DFSOS-1

DFSOS-2

APG

ADMM

SVM

1-NN

5-NN

10-NN

Arrowhead

Accuracy Run-time (s) Cardinality

0.701 (0.0) 0.924 (0.695) 0.554 (0.066)

0.701 (0.0) 0.819 (0.243) 0.554 (0.066)

0.688 (0.0) 0.67 (0.259) 0.891 (0.006)

0.71 (0.0) 0.271 (0.04) 0.572 (0.004)

0.754 (0.0) 0.362 (0.499) –

0.76 (0.0) 0.026 (0.003) –

0.714 (0.0) 0.017 (0.0) –

0.566 (0.0) 0.013 (0.0) –

Beef

Accuracy Run-time (s) Cardinality

0.867 (0.0) 0.908 (0.726) 0.979 (0.0)

0.867 (0.0) 0.842 (0.285) 0.979 (0.0)

0.579 (0.0) 1.343 (0.665) 0.918 (0.007)

0.642 (0.001) 0.459 (0.131) 0.587 (0.012)

0.867 (0.0) 0.364 (0.394) –

0.7 (0.0) 0.029 (0.002) –

0.633 (0.0) 0.019 (0.0) –

0.433 (0.0) 0.017 (0.0) –

BME

Accuracy Run-time (s) Cardinality

0.896 (0.0) 0.955 (0.692) 0.271 (0.029)

0.896 (0.0) 0.725 (0.19) 0.271 (0.029)

0.924 (0.0) 0.275 (0.061) 0.463 (0.025)

0.926 (0.0) 0.231 (0.029) 0.404 (0.004)

0.893 (0.0) 0.405 (0.645) –

0.84 (0.0) 0.04 (0.007) –

0.74 (0.0) 0.026 (0.002) –

0.573 (0.0) 0.015 (0.0) –

Car

Accuracy Run-time (s) Cardinality

0.77 (0.001) 1.018 (0.514) 0.793 (0.085)

0.77 (0.001) 0.937 (0.268) 0.793 (0.085)

0.75 (0.0) 0.971 (0.491) 0.977 (0.0)

0.765 (0.0) 0.48 (0.184) 0.612 (0.003)

0.817 (0.0) 0.338 (0.429) –

0.717 (0.0) 0.025 (0.003) –

0.583 (0.0) 0.023 (0.001) –

0.6 (0.0) 0.016 (0.001) –

CBF

Accuracy Run-time (s) Cardinality

0.688 (0.003) 0.713 (0.234) 0.294 (0.036)

0.688 (0.003) 0.669 (0.171) 0.294 (0.036)

0.852 (0.0) 0.498 (0.302) 0.282 (0.002)

0.851 (0.0) 0.285 (0.068) 0.309 (0.004)

0.862 (0.0) 0.442 (0.723) –

0.844 (0.0) 0.049 (0.015) –

0.764 (0.0) 0.032 (0.003) –

0.649 (0.0) 0.02 (0.001) –

EthanolLevel

Accuracy Run-time (s) Cardinality

0.74 (0.0) 2.098 (0.647) 0.997 (0.0)

0.74 (0.0) 1.946 (0.307) 0.997 (0.0)

0.391 (0.001) 13.223 (0.826) 0.952 (0.0)

0.417 (0.001) 6.632 (0.847) 0.827 (0.005)

0.838 (0.0) 1.284 (0.105) –

0.334 (0.0) 0.064 (0.002) –

0.356 (0.0) 0.061 (0.001) –

0.334 (0.0) 0.06 (0.001) –

Fish

Accuracy Run-time (s) Cardinality

0.823 (0.0) 0.533 (0.096) 0.997 (0.0)

0.823 (0.0) 0.571 (0.14) 0.997 (0.0)

0.802 (0.0) 3.696 (2.815) 0.97 (0.0)

0.797 (0.0) 1.425 (0.256) 0.759 (0.001)

0.834 (0.0) 0.413 (0.108) –

0.811 (0.0) 0.03 (0.002) –

0.766 (0.0) 0.025 (0.0) –

0.754 (0.0) 0.018 (0.0) –

Mallat

Accuracy Run-time (s) Cardinality

0.787 (0.0) 14.538 (1.918) 0.963 (0.0)

0.709 (0.001) 13.451 (0.356) 0.719 (0.001)

0.919 (0.0) 16.903 (1.541) 0.918 (0.0)

0.93 (0.0) 9.987 (1.241) 0.463 (0.004)

0.888 (0.0) 0.491 (0.103) –

0.738 (0.0) 0.027 (0.002) –

0.623 (0.0) 0.022 (0.001) –

0.5 (0.0) 0.016 (0.001) –

Meat

Accuracy Run-time (s) Cardinality

0.95 (0.0) 0.508 (0.331) 0.997 (0.0)

0.95 (0.0) 0.444 (0.12) 0.997 (0.0)

0.983 (0.0) 0.829 (0.501) 0.956 (0.0)

0.983 (0.0) 0.366 (0.079) 0.523 (0.008)

0.933 (0.0) 0.395 (0.586) –

0.933 (0.0) 0.05 (0.016) –

0.883 (0.0) 0.023 (0.001) –

0.883 (0.0) 0.012 (0.0) –

OliveOil

Accuracy Run-time (s) Cardinality

0.867 (0.0) 0.202 (0.012) 0.995 (0.0)

0.867 (0.0) 0.241 (0.023) 0.995 (0.0)

0.9 (0.0) 0.922 (0.546) 0.944 (0.0)

0.9 (0.0) 0.505 (0.2) 0.596 (0.005)

0.867 (0.0) 0.547 (1.269) –

0.833 (0.0) 0.048 (0.01) –

0.8 (0.0) 0.034 (0.004) –

0.767 (0.0) 0.015 (0.0) –

Plane

Accuracy Run-time (s) Cardinality

0.978 (0.0) 1.624 (1.753) 0.518 (0.026)

0.978 (0.0) 1.282 (0.413) 0.518 (0.026)

0.981 (0.0) 3.093 (1.615) 0.623 (0.003)

0.981 (0.0) 0.731 (0.137) 0.573 (0.003)

0.981 (0.0) 0.462 (0.279) –

0.962 (0.0) 0.028 (0.003) –

0.952 (0.0) 0.022 (0.0) –

0.952 (0.0) 0.016 (0.0) –

Rock

Accuracy Run-time (s) Cardinality

0.66 (0.0) 2.255 (1.919) 0.199 (0.001)

0.66 (0.0) 1.961 (0.178) 0.199 (0.001)

0.673 (0.019) 2.911 (0.565) 0.552 (0.158)

0.782 (0.007) 1.858 (0.075) 0.064 (0.001)

0.84 (0.0) 0.251 (0.134) –

0.76 (0.0) 0.03 (0.003) –

0.42 (0.0) 0.021 (0.001) –

0.46 (0.0) 0.016 (0.0) –

StarLightCurves

Accuracy Run-time (s) Cardinality

0.774 (0.0) 18.869 (19.86) 0.478 (0.086)

0.774 (0.0) 17.648 (12.908) 0.478 (0.086)

0.774 (0.0) 4.218 (0.281) 0.822 (0.003)

0.769 (0.0) 11.876 (6.043) 0.399 (0.008)

0.926 (0.0) 3.034 (0.102) –

0.844 (0.0) 0.074 (0.003) –

0.837 (0.0) 0.068 (0.001) –

0.816 (0.0) 0.06 (0.0) –

Symbols

Accuracy Run-time (s) Cardinality

0.824 (0.0) 4.107 (3.751) 0.179 (0.0)

0.827 (0.001) 3.287 (0.892) 0.206 (0.0)

0.796 (0.002) 3.296 (0.527) 0.783 (0.069)

0.824 (0.001) 1.779 (0.34) 0.463 (0.023)

0.837 (0.0) 0.319 (0.089) –

0.885 (0.0) 0.023 (0.002) –

0.677 (0.0) 0.019 (0.001) –

0.183 (0.0) 0.014 (0.001) –

UWaveGestureLibraryALL

Accuracy Run-time (s) Cardinality

0.5 (0.0) 1.571 (0.257) 0.998 (0.0)

0.5 (0.0) 1.282 (0.168) 0.998 (0.0)

0.834 (0.0) 20.479 (0.884) 0.946 (0.0)

0.828 (0.0) 23.313 (1.142) 0.736 (0.0)

0.881 (0.0) 1.034 (0.122) –

0.949 (0.0) 0.059 (0.002) –

0.94 (0.0) 0.055 (0.001) –

0.931 (0.0) 0.048 (0.001) –

Table 2: Out-of-sample prediction accuracy, run-time in seconds, and fraction of non-zero features for discriminant vectors and classifiers for second subset of data from the UCR Time Series repository. 4.3.2

Classification of UCR Time-Series Data

Summary statistics for our analysis of the UCR time-series data can be found in Table 2. As before, we performed one-sided Wilcoxon tests to compare pairs of classification methods in terms of classification error, run-time, and cardinality. Figure 7 summarizes the results of these hypothesis tests. We observe similar outcomes to that of the experiments involving Gaussian data. DFSOS has similar run-time and classification performance to that of the two ASDA methods, especially ADMM. However, inspecting Table 2 suggests that when DFSOS does provide an improvement in classification accuracy over ASDA then this improvement can be substantial, e.g., for the Beef and EthanolLevel data. In this case, we can infer that the improved stability obtained by eliminating deflationary errors allows computation of more accurate discriminant-based classifiers.

Conclusion

We have proposed an algorithmic framework for deflation-free sparse optimal scoring (DFSOS), a novel approach to sparse discriminant analysis that addresses fundamental limitations of existing deflationbased methods. By reformulating the sparse optimal scoring problem with an explicit global orthogonality constraint, DFSOS computes all discriminant vectors simultaneously rather than sequentially, eliminating error propagation inherent in deflation-based approaches. The proposed algorithm decomposes the problem into tractable subproblems for scoring vectors, dis-

$'00

690

$3 *

69 0

') 62 6

(a) Classification Accuracy

(b) Run-times

$3 *

') 62 6

')626

$3*

$'00

') 62 6

$3 *

') 62 6

1XPEHURIH0UHMHFWLRQV:LOFR[RQ&DUGLQDOLW\

')626

$3*

$'00

') 62 6

1XPEHURIH0UHMHFWLRQV:LOFR[RQ5XQ7LPH

')626 ')626

1XPEHURIH0UHMHFWLRQV:LOFR[RQ$FFXUDF\

Figure 7: Number of time-series datasets from UCR repository for which we observe statistically significant improvement in classification accuracy, computational efficiency, and cardinality of discriminants when using row method compared to column method. That is, the (i, j)-entry is the number of data sets where we reject the null hypothesis H0 : there is no difference between method i and method j using the Wilcoxon test with one-sided alternative hypothesis Ha : method i is better than method j.

(l) Rock

(m) StarLightCurves

1 1

(j) OliveOil $YHUDJH&RVLQH6LPLODULW\8:DYH*HVWXUH/LEUDU\$OO

(n) Symbols

69 0

$' 00

69 0

$3 *

$' 00

62 6

')626 ')626

$3*

$'00

690

62 6

69 0

$' 00

69 0

$' 00

$3 *

62 6

')626 ')626

690

$3*

$'00

690

$YHUDJH&RVLQH6LPLODULW\6\PEROV

(i) Meat

$'00

69 0

$' 00

$3*

$3 *

62 6

')626 ')626

$3*

$'00

690

$YHUDJH&RVLQH6LPLODULW\6WDU/LJKW&XUYHV

69 0

$' 00

$3 *

(h) Mallat

62 6

$YHUDJH&RVLQH6LPLODULW\5RFN

(g) Fish

$YHUDJH&RVLQH6LPLODULW\3ODQH

(k) Plane

$3 *

62 6

(f) EthanolLevel

62 6

')626 ')626

69 0

$3*

$' 00

$'00

$3 *

690

62 6

')626 ')626

$3*

$'00

690

62 6

$YHUDJH&RVLQH6LPLODULW\2OLYH2LO

(e) CBF

69 0

$' 00

$3 *

62 6

')626 ')626

690

$3*

$'00

690

$'00

69 0

$' 00

$3*

$3 *

62 6

')626 ')626

$3*

$'00

690

$YHUDJH&RVLQH6LPLODULW\0HDW

(d) Car

$YHUDJH&RVLQH6LPLODULW\0DOODW

69 0

$3 *

$' 00

') 6

62 6

')626 ')626

$3*

$'00

690

') 6

')626 ')626

$3*

$'00

690

62 6

1 1

26 ') 6

')626 ')626

$YHUDJH&RVLQH6LPLODULW\)LVK

62 6

1 1

26 ') 6

(b) Beef

')626 ')626

$YHUDJH&RVLQH6LPLODULW\(WKDQRO/HYHO

62 6

(a) Arrowhead

') 6

')626 ')626

$3*

$'00

690

') 6

$YHUDJH&RVLQH6LPLODULW\&%)

1 1

') 6

')626 ')626

$3*

$'00

690

') 6

$YHUDJH&RVLQH6LPLODULW\&DU

1 1

$'00

1 1

$3*

') 6

')626 ')626

$3*

$'00

690

$YHUDJH&RVLQH6LPLODULW\%0(

$YHUDJH&RVLQH6LPLODULW\%HHI

')626 ')626

$YHUDJH&RVLQH6LPLODULW\$UURZKHDG

(o) UWaveGestureLibraryAll

Figure 8: Average cosine similarity between out-of-sample predictions for each pair of methods. The diagonal entries indicate the average cosine similarity between all pairs of prediction vectors made by each method across all 10 repetitions for each data set.

criminant vectors, and orthogonality enforcement. We establish convergence of DFSOS to stationary points of the augmented Lagrangian under mild conditions on the penalty parameter, providing theoretical grounding for the method. Extensive numerical experiments on both synthetic Gaussian data and real-world time series from the UCR Time-Series Repository demonstrate that DFSOS achieves classification accuracy comparable to or better than existing deflation-based methods (e.g., APG and ADMM variants of ASDA), particularly in challenging settings with high feature correlation. While DFSOS requires modestly increased computational time compared to ASDA methods (typically 1-2 seconds additional per trial), it delivers statistically significant improvements in classification accuracy across most experimental settings. The method produces slightly denser discriminant vectors than ASDA, reflecting a favorable trade-off between sparsity and discriminative power. Notably, DFSOS achieves classification performance matching state-of-the-art support vector machines on linearly separable data while maintaining the interpretability advantages of sparse linear methods. The observed increase in computation used by DFSOS compared to the existing ASDA methods suggests an important avenue for potential research: the development of faster algorithms for solution of (8). The past few years have seen the development many alternatives to the splitting method of Lai and Osher [20] for orthogonality constrained optimization; see the recent manuscripts [1, 5, 29] for comprehensive treatments of optimization on manifolds. Specialization of these methods to the deflationfree sparse optimal scoring problem (8) could lead to improved computational efficiency. We plan a rigorous empirical comparison of specialization of a variety of state of the art algorithms for manifold optimization to deflation-free sparse optimal scoring. Moreover, further study of the convergence and consistency properties of the proposed DFSOS algorithm is needed; for example, analysis of the rate of convergence of DFSOS would greatly improve our understanding of the algorithm. The consistent improvements observed in our experiments, particularly for highly correlated data, suggest that deflation-free methods offer a more principled and effective approach to sparse discriminant analysis. By directly enforcing global orthogonality constraints rather than sequential deflation, DFSOS provides a coherent framework for feature selection and dimension reduction in high-dimensional classification problems. This demonstrates that the additional computational structure required to handle orthogonality constraints globally yields tangible benefits in classification performance and numerical stability, making deflation-free approaches a valuable addition to the sparse discriminant analysis toolbox.

Acknowledgments The authors acknowledge the use of the IRIDIS High Performance Computing Facility, and associated support services at the University of Southampton, in the completion of this work. B. Ames was supported in part by the National Science Foundation Grants #20212554 and #2108645, as well as the University of Alabama Research Grants RG14678 and RG14838.

References [1] P.-A. Absil, R. Mahony, and R. Sepulchre. Optimization algorithms on matrix manifolds. Princeton University Press, Princeton, NJ, 2008. [2] S. Atkins, G. Einarsson, L. Clemmensen, and B. Ames. Proximal methods for sparse optimal scoring and discriminant analysis. Advances in Data Analysis and Classification, 17(4):983–1036, 2023. [3] G. Bécigneul and O.-E. Ganea. arXiv:1810.00760, 2018.

Riemannian adaptive optimization methods.

arXiv preprint

[4] S. Bonnabel. Stochastic gradient descent on Riemannian manifolds. IEEE Transactions on Automatic Control, 58(9):2217–2229, 2013. [5] N. Boumal. An introduction to optimization on smooth manifolds. Cambridge University Press, Cambridge, 2023. [6] S. Burer and R. D. Monteiro. A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization. Mathematical programming, 95(2):329–357, 2003. [7] T. Cai and W. Liu. A direct estimation approach to sparse linear discriminant analysis. Journal of the American statistical association, 106(496):1566–1577, 2011. 19

[8] S. Chen, S. Ma, A. Man-Cho So, and T. Zhang. Nonsmooth optimization over the Stiefel manifold and beyond: Proximal gradient method and recent variants. SIAM Review, 66(2):319–352, 2024. [9] L. Clemmensen, T. Hastie, D. Witten, and B. Ersbøll. Sparse discriminant analysis. Technometrics, 53(4):406–413, 2011. [10] H. A. Dau, A. Bagnall, K. Kamgar, C.-C. M. Yeh, Y. Zhu, S. Gharghabi, C. A. Ratanamahatana, and E. Keogh. The UCR time series archive. IEEE/CAA Journal of Automatica Sinica, 6(6):1293– 1305, 2019. [11] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani. Least angle regression. Annals of Statistics, pages 407–451, 2004. [12] J. Fan, Y. Feng, and X. Tong. A road to classification in high dimensional space: the regularized optimal affine discriminant. Journal of the Royal Statistical Society Series B: Statistical Methodology, 74(4):745–771, 2012. [13] B. Gao, N. T. Son, P.-A. Absil, and T. Stykel. Riemannian optimization on the symplectic Stiefel manifold. SIAM Journal on Optimization, 31(2):1546–1575, 2021. [14] Y. Guo, T. Hastie, and R. Tibshirani. Regularized linear discriminant analysis and its application in microarrays. Biostatistics, 8(1):86–100, 2007. [15] T. Hastie, A. Buja, and R. Tibshirani. Penalized discriminant analysis. The Annals of Statistics, 23(1):73–102, 1995. [16] T. Hastie, R. Tibshirani, and A. Buja. Flexible discriminant analysis by optimal scoring. Journal of the American statistical association, 89(428):1255–1270, 1994. [17] T. Hastie, R. Tibshirani, and J. H. Friedman. The elements of statistical learning: data mining, inference, and prediction, volume 2. Springer, 2009. [18] C. He, J. Li, B. Jiang, S. Ma, and S. Zhang. On relatively smooth optimization over Riemannian manifolds. arXiv preprint arXiv:2508.03048, 2025. [19] B. Jiang, X. Meng, Z. Wen, and X. Chen. An exact penalty approach for optimization with nonnegative orthogonality constraints. Mathematical Programming, 198(1):855–897, 2023. [20] R. Lai and S. Osher. A splitting method for orthogonality constrained problems. Journal of Scientific Computing, 58:431–449, 2014. [21] Z. Lai, L.-H. Lim, and T. Tang. Stiefel optimization is NP-hard. arXiv preprint arXiv:2507.02839, 2025. [22] Z. Lai, L.-H. Lim, and K. Ye. Grassmannian optimization is NP-hard. SIAM Journal on Optimization, 35(3):1939–1962, 2025. [23] X. Li, S. Chen, Z. Deng, Q. Qu, Z. Zhu, and A. Man-Cho So. Weakly convex optimization over Stiefel manifold using Riemannian subgradient-type methods. SIAM Journal on Optimization, 31(3):1605–1634, 2021. [24] Z. Lin, H. Li, and C. Fang. Alternating direction method of multipliers for machine learning. Springer, Singapore, 2022. [25] C. Liu and N. Boumal. Simple algorithms for optimization on Riemannian manifolds with constraints. Applied Mathematics & Optimization, 82(3):949–981, 2020. [26] Q. Mai and H. Zou. A note on the connection and equivalence of three sparse linear discriminant analysis methods. Technometrics, 55(2):243–246, 2013. [27] Q. Mai, H. Zou, and M. Yuan. A direct approach to sparse discriminant analysis in ultra-high dimensions. Biometrika, 99(1):29–42, 2012. [28] S. Rosset and J. Zhu. Piecewise linear regularized solution paths. The Annals of Statistics, pages 1012–1030, 2007. 20

[29] H. Sato. Riemannian optimization and its applications, volume 670. Springer, Berlin, 2021. [30] J. Shao, Y. Wang, X. Deng, and S. Wang. Sparse linear discriminant analysis by thresholding for high dimensional data. The Annals of Statistics, 39(2):1241–1265, 2011. [31] The MathWorks Inc. Statistics and Machine Learning Toolbox, 2022. [32] The MathWorks Inc. MATLAB version: 9.14.0 (R2023a), 2023. [33] The MathWorks Inc. Statistics and Machine Learning Toolbox (R2023a), 2023. [34] Y. Wang, W. Yin, and J. Zeng. Global convergence of admm in nonconvex nonsmooth optimization. Journal of Scientific Computing, 78:29–63, 2019. [35] Z. Wen and W. Yin. A feasible method for optimization with orthogonality constraints. Mathematical Programming, 142(1):397–434, 2013. [36] D. M. Witten and R. Tibshirani. Penalized classification using Fisher’s linear discriminant. Journal of the Royal Statistical Society Series B: Statistical Methodology, 73(5):753–772, 2011. [37] H. Zhang and S. Sra. First-order methods for geodesically convex optimization. In Conference on learning theory, pages 1617–1638. PMLR, 2016. [38] H. Zou and T. Hastie. Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society Series B: Statistical Methodology, 67(2):301–320, 2005.

Record · ID 141484 · SHA-256 d3be11c68c2541a7
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.