Interpretable nonconvex submodule clustering algorithm using ℓr-induced tensor nuclear norm and ℓ2,p column sparse norm with global convergence guarantees

TL;DR

This paper introduces a new clustering method for 2D images that preserves data structure and improves performance.

Contribution

The novel 2D-NLRSC method uses tensor nuclear norm and ℓ2,p regularization for nonconvex clustering with global convergence guarantees.

Findings

2D-NLRSC outperforms existing methods on real image datasets.

The ℓr-induced tensor nuclear norm accurately approximates tensor rank.

The framework preserves higher-order correlations in 2D image data.

Abstract

Tensor-based subspace clustering algorithms have garnered significant attention for their high efficiency in clustering high-dimensional data. However, when dealing with 2D image data, traditional vectorization operations in most algorithms tend to undermine the correlations of higher-order tensor terms. To tackle this limitation, this paper proposes a non-convex submodule clustering approach (2D-NLRSC) that leverages sparse and low-rank representations for 2D image data. An ℓr-induced tensor nuclear norm is introduced to approximate the tensor rank precisely. Instead of vectorizing each 2D image, the framework arranges samples as lateral slices of a third-order tensor. It employs the t-product operation to generate an optimal representation tensor with low-rank constraint. The proposed method combines ℓq-norm induced clustering awareness with laplacian regularization to obtain a representation tensor with a diagonal structure. Additionally, 2D-NLRSC incorporates the ℓ2,p-norm as a regularization term, taking advantage of its excellent invariance, continuity, and differentiability. Experimental results on real image datasets validate the superior performance of the 2D-NLRSC model.

Figures

The framework of proposed 2D-NLRSC.

Comparison of ACC and NMI for Different Algorithms On COIL-20 Dataset.

Ablation study of the proposed method: Comparisons on dataset (A) CMU-PIE and dataset (B) Yale.

The performance of 2D-NLRSC on the ORL and Yale datasets under different K values.

ACC and NMI of the 2D-NLRSC method on ORL and Yale datasets: Dependence on parameters r,p,q and λ2 ((A),(E): ORL-r,p; (B),(F): Yale-r,p; (C),(G): ORL-q,λ2; (D),(H): Yale-q,λ2).

ACC and NMI of the 2D-NLRSC method on JAFFE, CMU-PIE, Yale and MNIST (L = 10) datasets: Dependence on parameters λ1 and λ3 ((A),(E): JAFFE; (B),(F): CMU-PIE; (C),(G): Yale; (D),(H): MNIST (L = 10)).

Error convergence curves of 2D-NLRSC on different datasets ((A): ORL; (B): JAFFE; (C): CMU-PIE; (D): Yale; (E): MNIST (L = 3); (F): MNIST (L = 5)).

Equations

Full text

1 Introduction

High-dimensional data often resides within lower-dimensional subspaces, revealing underlying subspace structures [1]. This insight underpins subspace clustering, which assumes that data points are sampled from multiple subspaces, with each point being linearly representable by a subset of points within the same subspace. The objective of subspace clustering is to divide data into clusters, where the points within each cluster lie in the same subspace. To address these challenges, a range of subspace clustering techniques have been proposed. Traditional subspace clustering methods, for example, Sparse Subspace Clustering (SSC) [2] and Low-Rank Representation (LRR) [3], have been extensively employed to partition data into multiple, potentially overlapping linear subspaces, thereby minimizing redundancy. These techniques have proven highly effective across various applications involving high-dimensional data, such as image clustering and temporal video segmentation [3]. Among the most prominent techniques, spectral clustering-based approaches have gained substantial traction over the past decade. These approaches involve constructing an affinity matrix and applying algorithms such as K-means or Normalized Cuts (Ncut). The versatility of these algorithms in diverse clustering tasks has been a key factor in their growing popularity. Recently, several works have focused on addressing performance degradation caused by noise in multi-view clustering. For example, Zhou et al. [4] proposed the RTOSNMF method, which tackles the issues of noise interference and insufficient utilization of inter-layer relationships in multi-layer network community detection by performing denoising via linear separation and the $[eqn]$ -norm, as well as exploring inter-layer relationships through nuclear norm-constrained low-rank property. Che et al. [5] put forward a robust multi-view clustering method based on weighted low-rank tensor approximation and noise separation: on one hand, it designs differentiated constraints for different types of noise to achieve fine grained noise elimination; on the other hand, it efficiently explores high-order correlations among multiple views through weighted low-rank tensor modeling. Xie et al. [6] reformulated the min-cut problem as a bi-bounded constraint problem and developed an algorithm suitable for size-constrained min-cut scenarios and generalizable to broader bi-bounded nonlinear optimal transport problems. Shi et al. [7] proposed an optimal transport framework with both upper and lower bound constraints to enhance clustering and classification by better capturing structural relationships in data.

However, conventional subspace clustering methods encounter challenges when dealing with high-dimensional tensor data [2,3]. This challenge lies at the core of the current research. These techniques often lack robust theoretical guarantees for clustering accuracy, particularly due to the high dimensionality of the data, which requires substantial time and memory resources. Our framework addresses this gap by providing a robust, tailored solution for image data clustering.

When handling imaging data, traditional subspace clustering methods often vectorize data samples, enabling the application of algorithms designed for vectorial data. While effective in some cases, this vectorization process disrupts the intrinsic multidimensional structure of images, thereby reducing the reliability of subsequent analysis. It increases the dimensionality of the data, exacerbating the curse of dimensionality, and it also destroys the inherent high-dimensional structure of the data. To address this limitation, various works have leveraged multilinear algebra tools to preserve and exploit the spatial characteristics of imaging data [8–10]. One significant innovation in this domain is the t-product, a matrix-analogous multiplication operation for third-order tensors that improves the exploitation of their internal structure. In this framework, the representation matrices from various views are arranged into slices of a three-dimensional tensor [11–13]. Utilizing operations from multilinear and abstract algebra enables the more effective exploration of third-order tensor properties. The introduction of the t-product simplifies tensor-tensor multiplication, enabling more efficient and accurate modeling of image data [14,15].

For instance, Wu [16] proposed a tensor-based submodule clustering method for 2D imaging data, which organized samples as lateral slices of third-order tensors via t-product, integrated low-rank constraints, manifold regularization, and a unified ADMM-spectral clustering framework with nonlinear extensions. Francis et al. [17] developed a robust unsupervised tensor subspace clustering method using $[eqn]$ -regularized tensor nuclear norm (TNN) to preserve geometric structures, enable slice-wise sparse-low-rank decomposition for noise removal. Madathil et al. [18] introduced a noise-robust tensor clustering approach via reweighted nuclear norms for enhanced low-rank representation, $[eqn]$ -norm structured sparsity, and explicit noise separation, achieving high accuracy under severe corruption. Francis et al. [19] proposed a single-stage framework for incomplete imaging data, unifying tensor clustering (via sparse t-linear combinations with mode-3 low-multirank constraints) and missing-data reconstruction (through low-rank lateral slice approximations), thus preserving spatial structure.

Inspired by the above, we adopt the t-product, based on circular convolution, to model the dynamic characteristics of consecutive image sequences. Using the t-product, image data samples are grouped into a third-order tensor and represented by our proposed tensor low-rank model, built around the union of free submodules. The resulting affinity information is then utilized for final clustering.

By combining t-product operations and tensor factorization, we extend the traditional LRR clustering method to accommodate multi-view data. This paper employs sparse and low-rank representations, similar to those mentioned above. However, it incorporates the latest developments in sparse coding to propose a highly interpretable $[eqn]$ -induced tensor nuclear norm. The optimization of this norm is based on the studies of non-convex $[eqn]$ -norm by Candes et al. [20], Zuo et al. [21], and Zha et al. [22], combining sparsity and low-rank properties.

Peng et al. [23] proposed a novel NMF algorithm with an $[eqn]$ -(pseudo) norm applied to the factor matrix to enforce sparse properties. Inspired by this novel column-sparse norm and influenced by [22], we naturally generalize the above definition along the frontal slice direction to the tensor-based $[eqn]$ -induced sparse approximation, thereby proposing an $[eqn]$ -(pseudo) norm.

A dissimilarity matrix M is constructed to impose constraints on the frontal slices of the representation tensor. Specifically, elements with smaller values (indicating higher similarity) in the dissimilarity matrix are used to inversely constrain the corresponding positions in the slices of the representation tensor to take on larger values. This approach enhances the representation tensor, yielding a more distinct cluster structure. We show that the nonconvex $[eqn]$ -norm, combined with our cluster-aware construction, effectively captures the block structure of the self-representation tensor $[eqn]$ . With an appropriate value of q, the $[eqn]$ -norm offers unique advantages in sparsity representation. In certain scenarios, selecting an appropriate value of q enables the $[eqn]$ -norm to achieve unique advantages in sparsity representation.

The key contributions of our paper are as follows.

The proposed 2D-NLRSC algorithm leverages the submodular self-expression property to address the issue in traditional tensor subspace clustering where vectorization of 2D images destroys higher-order tensor correlations.It directly constructs image samples as lateral slices of a third-order tensor, effectively avoiding structural information loss caused by vectorization.The $[eqn]$ -induced tensor nuclear norm is employed to accurately approximate the tensor rank, and combined with t-product operations to construct an optimal representation tensor under low-rank constraints, resolving the problem that traditional convex tensor nuclear norms fail to precisely characterize tensor low-rank structures.The $[eqn]$ -norm is incorporated into 2D-NLRSC as a regularization term, utilizing its excellent invariance, continuity, and differentiability to enhance the model’s resilience to outliers while optimizing the selection of representative data points to reduce redundancy. Additionally, the clustering-aware property of the $[eqn]$ -norm is combined with Laplacian regularization to guide the representation tensor in forming block-diagonal structures consistent with clustering objectives.An efficient alternation direction method of multipliers (ADMM) algorithm is proposed, with its convergence rigorously proven based on the Karush-Kuhn-Tucker (KKT) conditions. Experiments conducted on real-world datasets validate the effectiveness of the proposed method.

The remainder of this paper is organized as follows. Sect 2 introduces the preliminaries. Sect 3 reviews the related work. Sect 4 presents the proposed 2D-NLRSC method and its optimization procedure. Experimental results are reported in Section 5, while Sect 6 provides the convergence analysis. Finally, Sect 7 concludes the paper.

2 Notations and preliminary considerations

Here, the relevant definitions and preliminary notations for the variables are first presented. For the sake of conciseness, the related notations are summarized in Table 1.

Table 1: Notations.

Definition 1. (t-product [10]): Let $[eqn]$ and $[eqn]$ . Then the t-product $[eqn]$ is a tensor of size $[eqn]$ , i.e.,

$$[eqn]$$

where $[eqn]$ is defined as

$$[eqn]$$

The operator $[eqn]$ and its corresponding inverse operator $[eqn]$ are defined as

$$[eqn]$$

Theorem 1. (t-SVD [24]): For $[eqn]$ , the t-SVD of $[eqn]$ is given by $[eqn]$ where $[eqn]$ and $[eqn]$ are orthogonal tensors. $[eqn]$ is an f-diagonal tensor, and * denotes the t-product.

Definition 2. (tensor $[eqn]$ -norm [25]) For a matrix $[eqn]$ , the $[eqn]$ -norm is defined as

$$[eqn]$$

where $[eqn]$ . Particularly, considering p = 1, the $[eqn]$ -norm reduces to the $[eqn]$ -norm:

$$[eqn]$$

We then extend matrix $[eqn]$ -norm to tensor $[eqn]$ -norm as follows

$$[eqn]$$

3 Related work

Consider a data set $[eqn]$ that contains N images partitioned into m categories. Conventional subspace clustering pipelines vectorize each image $[eqn]$ into a flattened vector $[eqn]$ , constructing a data matrix $[eqn]$ . These vectors are hypothesized to reside near a union of m low-dimensional subspaces embedded in $[eqn]$ . The clustering objective is to group data points based on their intrinsic subspace affiliations.

Although the above vectorization approach has demonstrated excellent performance in numerous applications, it fails to account for the spatial structure of images. Additionally, these methods, which are typically used to represent the linear combinations of samples in subspaces, are unable to capture shifted copies within submodules. In contrast, the t-product [10] offers a novel algebraic approach that generalizes matrix multiplication to third-order tensors.

3.1 Linear algebra with the t-product

To preserve the spatial structure of data, N oriented matrices of size $[eqn]$ are stacked into a third-order tensor $[eqn]$ . Here, $[eqn]$ denotes the set of tube fibers of dimension $[eqn]$ , and $[eqn]$ denotes the set of oriented matrices of size $[eqn]$ . The goal is to define a multiplication operation between tube fibers, enabling “linear” combinations of oriented matrices in which the coefficients are themselves tube fibers rather than scalar values.

Following [26], the set $[eqn]$ can be regarded as a module over the ring $[eqn]$ . From this viewpoint, the t-product provides a natural generalization of matrix multiplication to third-order tensors, where the multiplication between elements is replaced by tube fiber multiplication, and the addition remains elementwise.

3.2 Representation of high-dimensional image data

Given an image matrix $[eqn]$ , we embed it into a third-order tensor by orienting it along the third mode

$$[eqn]$$

where each lateral slice $[eqn]$ is a tube fiber of length D. Let $[eqn]$ denote the collection of N such oriented tensors.

In the t-product framework, $[eqn]$ is viewed as a free module of rank n1 over the ring $[eqn]$ . Given a generating set (dictionary) $[eqn]$ , any $[eqn]$ admits a t-linear representation

$$[eqn]$$

For computational purposes, let $[eqn]$ and $[eqn]$ . Then the above expression can be compactly written as

$$[eqn]$$

Applying the discrete Fourier transform (DFT) along the third mode yields

$$[eqn]$$

where $[eqn]$ , $[eqn]$ , and $[eqn]$ denote the k-th frontal slices in the Fourier domain. This block-diagonal structure allows independent processing of each frequency slice, retaining spatial structure while enabling efficient computation.

By preserving the tensor structure of the original images, this representation mitigates the loss of spatial correlations caused by vectorization, and naturally supports linear modeling in the tensor algebra framework.

3.3 Submodule clustering by sparse and low-rank representation

Based on this theoretical framework, the sparse submodule clustering (SSmC) algorithm, as introduced in [26], can be expressed as the subsequent optimization problem

$$[eqn]$$

where $[eqn]$ denotes the representation tensor.

Inspired by the concept that images from different free submodules [27] should exhibit low correlation, Wu et al. proposed the structurally constrained low-rank submodular clustering (SCLRSmC) in [14]. The learning process for this method can be expressed as

$$[eqn]$$

where $[eqn]$ is a predefined weight matrix related to the data, and $[eqn]$ represents the Hadamard product.

Additionally, Wu [28] proposed an online low-rank tensor subspace clustering (OLRTSC) algorithm based on the nonconvex reformulation of tensor low-rank representation (TLRR) and the t-SVD framework, aiming to recover efficiently and cluster tensor data. This approach significantly reduces computational complexity and storage costs, handles dynamic data, and extends to scenarios with missing data. He also introduced an outlier-robust tensor low-rank representation (OR-TLRR) method, which simultaneously performs outlier detection and tensor data clustering within the t-SVD framework [29]. Research shows that t-SVD effectively reduces data dimensionality, extracts key features, and enhances data processing efficiency.

Although the above algorithms have shown promising results on real datasets, approximating tensor rank with the nuclear norm remains an imprecise problem. Since the nuclear norm treats every eigenvalue equally and penalizes noise components, it can lead to a suboptimal representation tensor. Therefore, this approach still requires improvement.

Recently, Jobin et al. [17,30] proposed a multi-view data clustering framework based on nonconvex low-rank tensor approximation

$$[eqn]$$

Definition 3. ( $[eqn]$ -induced tensor nuclear norm): Consider a tensor, $[eqn]$ with t-SVD $[eqn]$ , then its $[eqn]$ -induced tensor nuclear norm can be expressed as

$$[eqn]$$

Notably, this framework employs the $[eqn]$ -induced tensor nuclear norm (TNN) as a low tensor rank constraint, which enhances the low-rank property of the representation tensor. However, Zuo et al. [21] demonstrated in the sparse coding context that $[eqn]$ is not necessarily optimal, and a general $[eqn]$ should be considered to balance rank sparsity and numerical stability. The generic $[eqn]$ -norm minimization problem takes the form

$$[eqn]$$

where y is the input vector (e.g., a singular value vector), and $[eqn]$ with $[eqn]$ .

Problem (15) can be efficiently solved by the generalized soft-thresholding (GST) algorithm [21], which alternates between a gradient descent step and a generalized shrinkage step. Specifically, given y, the GST update is

$$[eqn]$$

where $[eqn]$ is applied element-wise:

$$[eqn]$$

with threshold $[eqn]$ and z^*^ being the positive root of $[eqn]$ .

As pointed out by Zha et al. [22], $[eqn]$ minimization with r<1 can more aggressively promote sparsity (and thus low rank) compared to the convex r = 1 case, often leading to better empirical performance in subspace clustering tasks.

Inspired by the above, we introduce the $[eqn]$ -induced tensor nuclear norm to enforce low-rank structure in the tensor. Unlike [30], our method allows adjusting the value of r between 0 and 1, enabling a more flexible and effective low-rank representation of the tensor.

Definition 4. ( $[eqn]$ -induced tensor nuclear norm): Consider a tensor $[eqn]$ with t-SVD $[eqn]$ . It’s $[eqn]$ -induced TNN can be defined as

$$[eqn]$$

where $[eqn]$ , $[eqn]$ and $[eqn]$ are the orthogonal tensors, and $[eqn]$ is an f-diagonal tensor whose frontal slices contain diagonal matrices, and $[eqn]$ denotes its Fourier transform.

Xie et al. proposed a nonconvex tensor multi-view clustering framework [25], which introduces a novel column-sparse norm instead of the squared Frobenius norm for $[eqn]$ : the $[eqn]$ - norm with $[eqn]$ . This norm exhibits properties such as invariance, continuity, and differentiability. Inspired by [25], this paper applies the $[eqn]$ - norm to 2D images.

4 Proposed method

4.1 Problem formulation and objective function

Recognizing that most conventional clustering methods rely on vectorized data representations and often overlook the intrinsic structure of image data, we focus on leveraging tensor representations, particularly for 2D images. Consequently, we propose a novel image submodule clustering method that preserves and leverages the inherent structure of image data to achieve more robust and reliable clustering outcomes.

For the 2D image samples X mentioned above, rather than vectorizing each image as in traditional subspace clustering methods, the images are arranged side by side to form a tensor $[eqn]$ . This tensor is then rotated to a dimension of $[eqn]$ , where $[eqn]$ denotes the size of the images and N indicates the number of images. Each image sample is represented using t-linear combinations. To obtain the optimal low-rank representation tensor $[eqn]$ , we apply an $[eqn]$ -induced tensor nuclear norm on $[eqn]$ . We add a sparse norm, the $[eqn]$ -norm, where $[eqn]$ , to strengthen the robustness of our approach. A slice-wise Laplacian regularization term centered around each image is also introduced to retain local details and expose nonlinear structures in high-dimensional space. This regularization guarantees consistent representations of data points within each image in the global space, enhancing overall efficiency. Therefore, the combination of sliced Laplacian, $[eqn]$ regularization, and tensor $[eqn]$ -induced TNN can be represented as follows

$$[eqn]$$

where $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ is the Laplacian matrix of $[eqn]$ , the adjacency matrix of a k-nearest neighbor (KNN) graph. It can be constructed by

$$[eqn]$$

Here, $[eqn]$ denotes the $[eqn]$ -neighborhood of $[eqn]$ under a given distance metric $[eqn]$ . Specifically, we define $[eqn]$ , where $[eqn]$ is derived by normalizing each lateral slice of $[eqn]$ such that $[eqn]$ for all $[eqn]$ .

The degree matrix $[eqn]$ is a diagonal matrix where the i-th diagonal element is computed as $[eqn]$ . Relevant details can be found in [25]. Furthermore, the suitable block diagonal structure representing tensors is beneficial for clustering multi-view data, enhancing the algorithm’s performance. In multi-view data, objects that belong to the same submodule exhibit a strong correlation, while objects from different submodules show relatively lower correlation [14]. To enforce a diagonal structure on $[eqn]$ , we prioritize images that exhibit a lower correlation between those belonging to different submodules. This correlation between different data points can be captured by the inconsistency weighted matrix $[eqn]$ . The elements of M are defined as

$$[eqn]$$

where $[eqn]$ represents the tensor that is acquired by normalizing every lateral slice of $[eqn]$ in a way that $[eqn]$ for $[eqn]$ , and $[eqn]$ is typically taken as the empirical average of all $[eqn]$ .

Unlike (12), we obtain a sparser solution by substituting the $[eqn]$ -norm with the $[eqn]$ -norm. By integrating all the above, our algorithm can be summarized as follows

$$[eqn]$$

where $[eqn]$ , $[eqn]$ , and $[eqn]$ are balance parameters. The flowchart of the 2D-NLRSC is shown in Fig 1. Our method processes the input image data as a third-order tensor, successfully preserving the original spatial structure of images and avoiding structural information loss due to vectorization. The first term employs the $[eqn]$ -induced tensor nuclear norm to impose a low-rank constraint on the representation tensor, thereby capturing the global subspace structure of the data. The second term applies the $[eqn]$ -norm to the Hadamard product of the representation tensor and the dissimilarity weighted matrix, promoting element-wise sparsity and forming a block-diagonal structure conducive to clustering. The third term uses the $[eqn]$ -norm to limit the error tensor, boosting column-wise sparsity and model robustness against outliers and noise. The fourth term uses Laplacian regularization to utilize the local geometric structure of the data, ensuring that the representation preserves neighborhood interactions within the intrinsic manifold and so strengthens intra-cluster cohesiveness and inter-cluster separation. .

The framework of proposed 2D-NLRSC.

By solving (22), the optimal representation tensor can be used to construct affinity matrices.The affinity matrix S is formulated as

$$[eqn]$$

Then, spectral clustering or the normalized graph cut algorithm is applied to the affinity matrix S to derive the final clustering outcomes.

4.2 Optimization

Three auxiliary variables $[eqn]$ , $[eqn]$ , $[eqn]$ are introduced into our optimization algorithm, which can be expressed as follows

$$[eqn]$$

The augmented Lagrangian function for this problem can be expressed as

$$[eqn]$$

1. $[eqn]$ Subproblem: Fixing all variables except $[eqn]$ , the update for $[eqn]$ is

$$[eqn]$$

Theorem 2. Let $[eqn]$ have the t-SVD decomposition $[eqn]$ . Consider the $[eqn]$ -induced tensor nuclear norm optimization problem

$$[eqn]$$

the optimal solution is given by

$$[eqn]$$

where $[eqn]$ is an f-diagonal tensor, whose diagonal elements are generated by the GST algorithm described in (16).

Applying Theorem 2, the optimal solution of (26) is

$$[eqn]$$

2. $[eqn]$ Subproblem: The update for $[eqn]$ is

$$[eqn]$$

Following the GST algorithm described in (16), the closed-form solution is

$$[eqn]$$

where $[eqn]$ .

3. $[eqn]$ Subproblem: The update for $[eqn]$ is

$$[eqn]$$

Setting the partial derivatives with respect to $[eqn]$ to zero gives

$$[eqn]$$

The optional solution for $[eqn]$ is derived as

$$[eqn]$$

4. $[eqn]$ Subproblem: The update rule for $[eqn]$ is formulated as

$$[eqn]$$

where $[eqn]$ , and $[eqn]$ is the v-th frontal slice of $[eqn]$ .

Let $[eqn]$ , $[eqn]$ be the symbol for the i-th column of E, and take $[eqn]$ to mean the i-th column of $[eqn]$ . The objective in (35) can be reformulated column-wise as

$$[eqn]$$

so that each $[eqn]$ can be solved separately. For a specific $[eqn]$ , the subproblem becomes

$$[eqn]$$

Here, $[eqn]$ can be treated as a special matrix, and a thin SVD can be applied. It follows that $[eqn]$ has exactly one singular value, given by

$$[eqn]$$

where $[eqn]$ represents the singular value of the input vector. Hence, the subproblem is equivalent to

$$[eqn]$$

Now, we introduce Lemma 1.

Lemma 1. [31] Consider two complex matrices A and $[eqn]$ . Let $[eqn]$ be defined as

$$[eqn]$$

where $[eqn]$ represents the vector consisting of the non-increasing singular values of A. If $[eqn]$ is a complex unitarily invariant function (or quasi-norm) and B has an SVD formulated as

$$[eqn]$$

then the optimal solution $[eqn]$ for the minimization problem

$$[eqn]$$

has the SVD $[eqn]$ , where

$$[eqn]$$

and

$$[eqn]$$

According to Lemma 1, the solution to (39) is expressed as

$$[eqn]$$

where ui and $[eqn]$ are respectively the left and right singular vectors corresponding to $[eqn]$ , and

$$[eqn]$$

which can be computed using the GST algorithm in (16).

It is evident that(refer to 7.2 of [32])

$$[eqn]$$

denotes a thin SVD of $[eqn]$ , where [1] represents a matrix with 1 as its only element. By replacing

$$[eqn]$$

$[eqn]$ into (45),we can derive $[eqn]$ .

5. $[eqn]$ Subproblem: The update for $[eqn]$ is

$$[eqn]$$

where $[eqn]$ . Using discrete Fourier transform, the solution is

$$[eqn]$$

Overall, the algorithm is summed up in Algorithm 1.

Algorithm 1 The algorithm of 2D-NLRSC.

Input: Given data tensor $[eqn]$ , dissimilarity matrix $[eqn]$ , and   parameters $[eqn]$ , $[eqn]$ , and $[eqn]$ .

Output: Representation tensor $[eqn]$ .

1: Initialization: $[eqn]$

$[eqn]$ , penalty parameter $[eqn]$ ,

$[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , and t = 0.

2: while not converge do

3:   Fix the others and update $[eqn]$ by (50);

4:   Fix the others and update $[eqn]$ by (45);

5:   Fix the others and update $[eqn]$ by (29);

6:   Fix the others and update $[eqn]$ by (31);

7:   Fix the others and update $[eqn]$ by (34);

8:   Fix the others and update the Lagrange multipliers $[eqn]$ , $[eqn]$ ,

$[eqn]$ and $[eqn]$ by

emsp;   $[eqn]$ ;

$[eqn]$ ;

$[eqn]$ ;

$[eqn]$ ;

9:   Update the parameter $[eqn]$ and $[eqn]$ by

$[eqn]$ ;

$[eqn]$ ;

10:   Check the convergence conditions if

$$[eqn]$$

11:   Output: Representation tensor $[eqn]$ .

12:   else, t = t + 1;

13: end while

4.3 Complexity analysis

The computational cost of our method primarily arises from the updates of $[eqn]$ and $[eqn]$ . When updating $[eqn]$ , it is necessary to compute the FFT and inverse FFT of an $[eqn]$ tensor along mode-3, as well as perform the SVD of $[eqn]$ matrices in the Fourier domain. These operations require $[eqn]$ computations per iteration. For updating $[eqn]$ , the process involves calculating the matrix inverse and the FFT and inverse FFT transformations. The computational complexity of these operations is $[eqn]$ . Consequently, the total computational complexity of the method is approximately $[eqn]$ , where T1 denotes the number of iterations required to solve (24) using ADMM.

5 Experiments

5.1 Datasets

Five image datasets are selected for this algorithm, and they are described in Table 2 as follows

1. ORL: The ORL dataset comprises 40 distinct subjects, each represented by 10 different images. The experimental settings align with those described in [16]. All 400 images are utilized, and each image is of a dimension of $[eqn]$ .2) JAFFE: The JAFFE dataset comprises 213 images of 7 facial expressions from ten Japanese female subjects. Ten participants were selected for this study, each contributing their first 20 images. All 200 images are resized to $[eqn]$ . The experimental setup follows the same approach outlined in [16].3) CMU-PIE: The CMU-PIE face dataset comprises 42,368 images of 68 subjects, exhibiting diverse poses, lighting conditions, and facial expressions. A subset of 735 images was created, containing the initial 49 images from each of the first 15 subjects. These images were subsequently resized to $[eqn]$ . The specific experiments and settings are in [16].4) Yale: The Yale face dataset consists of 165 images of 15 individuals, each with a size of $[eqn]$ . 11 distinct facial images, exhibiting varied expressions and lighting, were provided by each participant. These images were uniformly resized to $[eqn]$ . The experimental conditions are aligned with those described in [16].5) MNIST:The MNIST dataset contains 70,000 centered $[eqn]$ images of handwritten digits (0-9). A subset of 1,000 images was created by selecting the first 100 images of each digit. The specific experimental settings are described in [16]. For each value of $[eqn]$ , the clustering process was repeated 20 times with randomly selected categories. The average of these 20 clustering results was then used for evaluation.6) COIL-20:The COIL-20 dataset contains 1,440 images from 20 different object categories. Each object category has 72 images captured from different angles, and the backgrounds of the objects were removed during the shooting process. These images were processed and downsampled to a size of $[eqn]$ .

Table 2: Statistics of the image datasets in the experiments.

In our experiment, Accuracy (ACC) and Normalized Mutual Information (NMI) are utilized as performance evaluation metrics. Higher values for both ACC and NMI indicate superior clustering performance.

5.2 Compared clustering algorithms

To evaluate the performance of the algorithms proposed in this paper, we selected the state-of-the-art clustering algorithms as benchmarks for our experimental comparisons

LRR [3] is a clustering algorithm that reveals the latent structure in data by constructing a low-rank representation matrix.SSC [2] solves a sparse optimization program to derive the sparse representation of data points from other points.LSR [33] leverages data correlation to achieve subspace segmentation through a grouping effect that tends to cluster highly correlated data together, significantly improving segmentation accuracy while ensuring efficiency.SC-LRR [34] extends the standard LRR by introducing a predefined weight matrix to analyze the structure of multiple disjoint subspaces, breaking through the restriction of the standard LRR that requires subspaces to be independent.TSC [35] achieves clustering of noisy and incompletely observed high-dimensional data into a union of low-dimensional subspaces and outliers by thresholding the correlations between data points to obtain an adjacency matrix.S3C [36] learns both the affinity matrix and segmentation results simultaneously through a joint optimization framework, expressing each data point as a structured, sparse linear combination of other data points.KSSC [37] extends sparse subspace clustering to nonlinear manifolds via the kernel trick, enhancing clustering performance through nonlinear mappings.SSmC [26] integrates the t-product operation into the sparse subspace clustering framework, enhancing the model’s ability to capture local features of data through its convolutional structure.SCLRSmC [14] introduces a free submodules theory, enabling low-rank representation learning directly in the tensor space and avoiding the inherent loss of structural information associated with traditional vectorization preprocessing.CLLRSmC [16] considers the intrinsic manifold structure of data and integrates the two distinct stages of learning a low-rank representation tensor and performing spectral clustering into a unified optimization framework.KCLLRSmC [16] combines manifold regularization with kernel methods for manifold clustering, eliminating the need for explicit mapping of data into the feature space. It serves as a nonlinear extension of CLLRSmC.

5.3 Experiments results and analysis

The subsequent experiments use ACC and NMI as clustering evaluation metrics. Their precise definitions are provided in [16]. Each experiment is replicated 20 times, and the average results are recorded in Tables 3, 4, 5, and 6. The best results are highlighted in bold. As shown in the tables, the proposed model significantly outperforms other models regarding both ACC and NMI across the ORL, JAFFE, CMU-PIE, Yale, and MNIST datasets.

Table 3: Experimental results on ORL and JAFFE datasets.The parameters are set as λ1=0.001, λ2=0.001, λ3=0.1 on ORL; λ1=0.0001, λ2=0.001, λ3=0.0001 on JAFFE.

Table 4: Experimental results on CMU-PIE and Yale datasets.The parameters are set as λ1=0.0001, λ2=0.001, λ3=0.0001 on CMU-PIE; λ1=0.0001, λ2=0.001, λ3=0.0001 on Yale.

Table 5: Experimental results on MNIST dataset (L = 3 and L = 5).The parameters are set as λ1=0.01, λ2=0.001, λ3=10 on MNIST.

Table 6: Experimental results on MNIST dataset (L = 8 and L = 10).The parameters are set as λ1=0.01, λ2=0.001, λ3=10 on MNIST.

On the ORL dataset, our method outperforms all other models, achieving improvements of 0.50% in ACC and 1.50% in NMI compared to the suboptimal KCLLRSmC.

On the JAFFE dataset, as shown in Table 3, the proposed model achieves a perfect score of 100% for both ACC and NMI. In particular, our model outperforms the suboptimal KCLLRSmC model by 0.36% in ACC and 0.77% in NMI.

Experiments on the CMU-PIE dataset demonstrate that 2D-NLRSC achieves superior performance compared to existing methods, yielding the highest ACC and NMI scores. Specifically, it surpasses the suboptimal method by 2.32% in ACC and 0.54% in NMI.

On the Yale dataset, 2D-NLRSC achieved the highest ACC and NMI scores, surpassing the best-performing suboptimal method by 1.03% and 0.77%, respectively.

On the MNIST dataset, compared to all other methods, the proposed 2D-NLRSC shows significant improvements in ACC and NMI. When the number of categories L = 8, the NMI is slightly lower than those of KCLLRSmC and CLLRSmC, but outperforms other algorithms. Compared with KCLLRSmC, when L = 3, the ACC and NMI metrics are improved by 2.16% and 3.59%, respectively. Similarly, when L = 5 and L = 10, the ACC and NMI metrics are improved by {0.37%, 0.52%} and {0.08%, 0.03%}, respectively.

Experimental evaluations on the COIL-20 dataset, detailed in Fig 2, illustrate that the 2D-NLRSC method yields significant advancements in clustering performance, reflected by its elevated ACC and NMI scores.

Comparison of ACC and NMI for Different Algorithms On COIL-20 Dataset.

In summary, the proposed 2D-NLRSC method significantly improves almost all datasets. These results strongly validate that 2D-NLRSC can more effectively capture high-order correlation information in samples through tensor rank approximation and joint sparse regularization based on the $[eqn]$ -norm and the $[eqn]$ -norm.

5.4 Ablation studies

To validate the role of each regularization term in the model, ablation experiments were conducted. Specifically, partial regularization terms were removed from the model 2D-NLSRC, and parameters were adjusted to achieve the optimal performance of the model. On the CMU-PIE and Yale dataset, three groups of experiments were designed, with each experiment omitting one regularization term to derive three algorithms. Through parameter optimization for each algorithm, the impact of each regularization term on the model’s performance was evaluated. The specific experimental results are shown in Fig 3.

Ablation study of the proposed method: Comparisons on dataset (A) CMU-PIE and dataset (B) Yale.

According to Fig 3, the suggested 2D-NLSRC method outperforms all other ablation variants in terms of clustering performance, as measured by ACC and NMI. We observed similar performance deterioration across all incomplete configurations in systematic ablation trials with individual regularization factors removed. On the CMU-PIE dataset, the complete 2D-NLSRC model improves ACC by {0.86%, 15.42%, 0.64%, 18.19%} and NMI by {1.22%, 7.34%, 0.83%, 12.29%} over the ablation versions 2D-NLSRC( $[eqn]$ ), 2D-NLSRC( $[eqn]$ ), 2D-NLSRC( $[eqn]$ ), and 2D-NLSRC( $[eqn]$ ). The Yale dataset shows performance increases of {2.36%, 9.76%, 1.45%, 13.18%} in ACC and {1.84%, 8.59%, 0.61%, 9.02%} in NMI.

These empirical findings give significant support for the efficacy of multi-component collaboration in our approach. The performance loss found when $[eqn]$ emphasizes the necessity of the $[eqn]$ -norm regularization that works in conjunction with the inconsistency matrix M. This combination applies differentiated constraints to sample pairs from various submodules, promoting the construction of a distinct block-diagonal structure in the representation tensor. The performance drop at $[eqn]$ confirms that the $[eqn]$ -norm introduces crucial column-wise sparse constraints, effectively eliminating interference from redundant information and outliers. The performance degradation observed when $[eqn]$ demonstrates that the Laplacian regularization term effectively preserves the local geometric structure of the data, and its absence leads to disrupted neighborhood relationships that are crucial for maintaining clustering coherence The significant performance degradation in the dual-ablated scenario 2D-NLSRC ( $[eqn]$ ) highlights the complimentary nature of these two regularization techniques, which jointly contribute to the improved clustering accuracy of 2D-NLSRC on complex datasets.

5.5 Sensitivity analysis

To comprehensively evaluate the robustness and stability of the proposed algorithm, we conduct a sensitivity analysis focusing on two key aspects: the parameter K in K-Nearest Neighbors (KNN) and the initialization of core variables ( $[eqn]$ ).

5.5.1 Sensitivity analysis of KNN parameter K.

We set the range of K from 5 to 30 and conduct experiments on both ORL and Yale datasets, with the specific results shown in Fig 4. The experimental data indicates that: on the ORL dataset, when K is in the range [5, 30], the fluctuation range of the algorithm’s accuracy ACC is controlled within 1.5%, and the fluctuation of normalized mutual information NMI does not exceed 1%; on the Yale dataset, the ACC fluctuation range is within 4%, and the NMI fluctuation does not exceed 3.6% . This result fully demonstrates the algorithm’s robustness to changes in K . Further analysis reveals that the relationship between K and algorithm performance is not linear. Within the tested range, the experimental results indicate that setting K=5 yields superior performance. Therefore, K is set to 5 for the subsequent experiments.

The performance of 2D-NLRSC on the ORL and Yale datasets under different K values.

5.5.2 Sensitivity analysis of initialization of core variables (𝒵,𝒞,𝒬,𝒯).

To comprehensively investigate the impact of the initialization of core variables ( $[eqn]$ ) on the 2D-NLRSC algorithm, we design multiple sets of comparative experiments: we perform 5 repeated experiments for zero initialization and three types of random Gaussian initializations with different variances ( $[eqn]$ , $[eqn]$ , $[eqn]$ ), and the final results are summarized in Table 7 after taking the average.On the ORL dataset, both ACC and NMI of zero initialization are the highest among all initialization methods; on the Yale dataset, its ACC and NMI are also significantly better than various random initializations, indicating that zero initialization can provide a better initial starting point for the algorithm, enabling the model to converge to a solution of higher quality. As the variance of the Gaussian distribution increases from 0.1 to 1, the ACC on the ORL dataset decreases from 0.8974 to 0.8785, and the NMI decreases from 0.9474 to 0.9408; on the Yale dataset, the ACC decreases from 0.9318 to 0.9273, and the NMI decreases from 0.9408 to 0.9348. It can be seen that the performance of random initialization shows a downward trend as the variance increases, and the larger the variance, the more obvious the performance degradation, indicating that the algorithm has a certain sensitivity to the variance change of random initialization.

Table 7: Performance comparison of 2D-NLRSC under different initializations on ORL and Yale datasets.

For both the ORL dataset (with 41 iterations for all) and the Yale dataset (with 36 iterations for all), the number of convergence iterations of different initialization methods is completely consistent. This indicates that although initialization affects the quality of the solution, it has no significant interference with the convergence speed of the algorithm, reflecting the stability of the algorithm in terms of convergence efficiency. Therefore, based on the comprehensive experimental evidence that zero initialization consistently achieves superior performance on both datasets without impairing convergence speed, it is adopted as the initialization method of choice in this paper.

5.6 Parameter sensitivity

In our method, $[eqn]$ , $[eqn]$ , and $[eqn]$ serve as balancing parameters, whereas r, p, and q denote hyperparameters. We systematically investigated their impact on clustering performance (quantified by ACC and NMI) across five benchmark datasets. Through exhaustive grid search, optimal parameter configurations were empirically determined for each dataset, as presented in [Tables 3–6](#pone.0339534.t003 pone.0339534.t004 pone.0339534.t005 pone.0339534.t006).

For the ORL and Yale datasets, we fixed all parameters except r and p, which were varied within $[eqn]$ to isolate their effects. To analyze the joint influence of q and $[eqn]$ on Yale and ORL datasets, we set $[eqn]$ $[eqn]$ and q $[eqn]$ while keeping the other variables constant. For the JAFFE, CMU-PIE, Yale and MNIST (L = 10) datasets, $[eqn]$ is fixed at 0.001, and $[eqn]$ , $[eqn]$ were varied within $[eqn]$ to evaluate their regulatory effects. Visualizations of these parameter analyses are provided in Figs 5 and 6.

ACC and NMI of the 2D-NLRSC method on ORL and Yale datasets: Dependence on parameters r,p,q and λ2 ((A),(E): ORL-r,p; (B),(F): Yale-r,p; (C),(G): ORL-q,λ2; (D),(H): Yale-q,λ2).

ACC and NMI of the 2D-NLRSC method on JAFFE, CMU-PIE, Yale and MNIST (L = 10) datasets: Dependence on parameters λ1 and λ3 ((A),(E): JAFFE; (B),(F): CMU-PIE; (C),(G): Yale; (D),(H): MNIST (L = 10)).

As shown in Fig 5, the value of r does not always correlate positively with clustering performance. On the ORL and Yale datasets, when r approaches 1, the clustering effect significantly declines, whereas setting r around 0.5 yields the best results. For parameter p, a threshold effect is evident: when p<0.7, clustering performance is constrained. As p increases, the algorithm performance improves gradually, however, when p = 1, the $[eqn]$ -norm degenerates into the $[eqn]$ -norm, which weakens sparsity constraints and degrades ACC. Based on the experimental results, the optimal range for p is [0.7,1). Parameter q exhibits a relatively stable influence, with stable and superior clustering performance observed when its value is within [0.4,0.9]. To ensure generalization ability and clustering performance of the model across different datasets, this paper selects $[eqn]$ , $[eqn]$ , $[eqn]$ for the ORL, JAFFE, CMU-PIE, Yale, MNIST and COIL-20 datasets, respectively.

As shown in Fig 6, after parameter tuning, the ACC and NMI of our algorithm exhibit similar trends across different datasets. Notably, when $[eqn]$ , the clustering performance deteriorates significantly. Similarly, as shown in Fig 5, on the Yale dataset, the performance also degrades when $[eqn]$ . Based on empirical results, we recommend setting $[eqn]$ , $[eqn]$ , and $[eqn]$ (as $[eqn]$ has a relatively stable impact).

Additionally, the running times of the 2D-NLRSC on different datasets are provided. As clearly shown in the Table 8, our algorithm exhibits significantly shorter running times than the CLIRSmc and KCLIRSmc algorithms on the ORL, JAFFE, and CMU-PIE datasets. Furthermore, our algorithm achieves higher ACC and NMI values compared to the KCLIRSmc algorithm. Overall, our algorithm demonstrates an acceptable time complexity.

Table 8: Running time (sec) comparison on different datasets.

6 Convergence analysis

The convergence of the 2D-NLRSC is examined from two distinct perspectives in this section. First, a rigorous theoretical analysis is provided to establish the convergence properties of the proposed method. Second, the convergence behavior is demonstrated through plots of various variables, as shown in Fig 7.

Error convergence curves of 2D-NLRSC on different datasets ((A): ORL; (B): JAFFE; (C): CMU-PIE; (D): Yale; (E): MNIST (L = 3); (F): MNIST (L = 5)).

Lemma 2. [31] Suppose $[eqn]$ is portrayed as $[eqn]$ , where $[eqn]$ has the SVD $[eqn]$ , $[eqn]$ , and f is differentiable. The gradient of F(X) at X is

$$[eqn]$$

where $[eqn]$ .

Lemma 3. [38] If the matrix H is positive definite, then

$$[eqn]$$

where $[eqn]$ represents the inner product, and a stands for the smallest eigenvalue of H.

Lemma 4. [38] In the finite-dimensional Euclidean space, every bounded sequence of vectors has a subsequence that converges.

Lemma 5. [39] Consider the nonconvex optimization problem with $[eqn]$ regularization

$$[eqn]$$

where $[eqn]$ , $[eqn]$ , and $[eqn]$ . Assume the function f satisfies the following conditions: f has Lf-Lipschitz continuous gradient, i.e.,

$$[eqn]$$

and f is bounded below on $[eqn]$ .

Then the following properties hold

(i) If $[eqn]$ is a local minimizer of (53), then it satisfies the first-order stationary condition

$$[eqn]$$

where $[eqn]$ and $[eqn]$ .(ii) Let $[eqn]$ be a stationary point with $[eqn]$ for some initial point $[eqn]$ and small $[eqn]$ . Define $[eqn]$ and let $[eqn]$ . Then each nonzero component satisfies

$$[eqn]$$

Lemma 6. $[eqn]$ is bounded.

Proof: In order to achieve the minimization of $[eqn]$ at the $[eqn]$ -th step as illustrated in equation (26), the optimal solution $[eqn]$ must satisfy the following condition

$$[eqn]$$

where $[eqn]$ . The expression $[eqn]$ has a singularity near $[eqn]$ . To circumvent this, we propose an approximation for $[eqn]$ with $[eqn]$

$$[eqn]$$

Let $[eqn]$ be expressed as $[eqn]$ . From Definition 4 and Lemma 2, it can be deduced that

$$[eqn]$$

Then we have

$$[eqn]$$

Therefore, it can be seen that $[eqn]$ is bounded. We can denote $[eqn]$ , where $[eqn]$ is the DFT matrix of size $[eqn]$ , and $[eqn]$ is its conjugate transpose. Given that $[eqn]$ , the following result is derived using the chain rule of matrix calculus

$$[eqn]$$

is bounded.

From the relations $[eqn]$ and $[eqn]$ , it follows that $[eqn]$ is bounded. $[eqn]$

Lemma 7. $[eqn]$ is bounded.

Proof: To prove that $[eqn]$ is bounded, consider the minimization of $[eqn]$ at the $[eqn]$ -th step as illustrated in (30). The optimal solution $[eqn]$ must satisfy

$$[eqn]$$

From the update $[eqn]$ , it follows that

$$[eqn]$$

The boundedness of the $[eqn]$ norm with power q is bounded [40] implies the boundedness of $[eqn]$ . $[eqn]$

Lemma 8. $[eqn]$ is bounded.

Proof: To achieve the minimization of $[eqn]$ at the $[eqn]$ -th step as illustrated in (32), the optimal solution $[eqn]$ requires to satisfy the following

$$[eqn]$$

due to the $[eqn]$ update by $[eqn]$ , then we can obtain

$$[eqn]$$

Based on the proof above, auxiliary variable $[eqn]$ is bounded, which implies that $[eqn]$ and its associated auxiliary variable $[eqn]$ are also bounded. According to the triangle inequality, it can be obtained that

$$[eqn]$$

Clearly, $[eqn]$ is bounded. $[eqn]$

Lemma 9. $[eqn]$ is bounded.

Proof: Consider the =of $[eqn]$ at the $[eqn]$ -th step as illustrated in (35). The optimal solution $[eqn]$ must satisfy

$$[eqn]$$

Given $[eqn]$ , the partial derivative of $[eqn]$ with respect to $[eqn]$ is expressed as follows

$$[eqn]$$

Thus, the following inequality holds

$$[eqn]$$

From the update $[eqn]$ , this can be rewritten as

$$[eqn]$$

Thus, $[eqn]$ is bounded. $[eqn]$

Theorem 3. Let $[eqn]$ be the sequence generated by Algorithm 1. Then the sequence $[eqn]$ satisfies the following two principles:

$[eqn]$ is bounded.Any accumulation point of $[eqn]$ is a KKT point of (25).

Proof: 1) Proof of the first part of Theorem 3:

Given the following update rules

$$[eqn]$$

It can be deduced that

$$[eqn]$$

Summing both sides of Eq (72) from t = 1 to k, we have

$$[eqn]$$

Since the sequences $[eqn]$ and $[eqn]$ are bounded, and

$$[eqn]$$

It can be found that the right of (73) is finite. Thus $[eqn]$ is bounded. We observed

$$[eqn]$$

According to Lemma 2, we can find that the right-hand side of (75) is bounded. Since the right of (75) is nonnegative, every term is also bounded. The boundedness of $[eqn]$ means that all singular values of $[eqn]$ are bounded. Consequently, this ensures that $[eqn]$ (the sum of the squares of the singular values) is also bounded. Therefore, the sequence $[eqn]$ is bounded.

According to the Lemma 3,

$$[eqn]$$

where $[eqn]$ is the smallest positive eigenvalue of the positive definite matrix L. Considering the expression $[eqn]$ is bounded, we can further deduce that $[eqn]$ is bounded.

Given that $[eqn]$ for all $[eqn]$ with constant M>0 and $[eqn]$ , we prove the sequence $[eqn]$ is bounded. For $[eqn]$ , the $[eqn]$ - norm is defined as

$$[eqn]$$

where $[eqn]$ satisfies

$$[eqn]$$

Define $[eqn]$ for each column j. Then $[eqn]$ . As $[eqn]$ and p > 0, for all j,

$$[eqn]$$

Since $[eqn]$ preserves element values, $[eqn]$ corresponds to the sum of squares over all entries of $[eqn]$ with fixed second index i2 = j,

$$[eqn]$$

Thus, for any entry $[eqn]$ ,

$$[eqn]$$

Set $[eqn]$ . Then $[eqn]$ for all $[eqn]$ . Hence, $[eqn]$ is bounded.

Similarly, the boundedness of $[eqn]$ implies that the sequences $[eqn]$ is also bounded. We have

$$[eqn]$$

because $[eqn]$ and $[eqn]$ are bounded, we deduce that $[eqn]$ is also bounded.

In conclusion, it can be proven that $[eqn]$ is bounded.

2) Proof of the second part of Theorem 3:

Above we have proved that the sequence $[eqn]$ generated by Algorithm 1 is ensured to be bounded. According to Lemma 4, which asserts that any bounded sequence in $[eqn]$ possesses a convergent subsequence, it follows that $[eqn]$ is bound to have at least a single point of accumulation. Let one of these points be denoted as as $[eqn]$ . Suppose, for the sake of generality, that $[eqn]$ .

Based on the update rule regarding $[eqn]$ , it can be concluded that $[eqn]$ . Then, by taking the limit of both sides of this equation

$$[eqn]$$

we obtain $[eqn]$ .

Similarly, we obtain $[eqn]$ . Based on the update rule for $[eqn]$ , it is evident that

$$[eqn]$$

then we can obtain $[eqn]$

Considering the $[eqn]$ -subproblem, it is observed that

$$[eqn]$$

Considering the $[eqn]$ -subproblem, it can be obtained that

$$[eqn]$$

Considering the $[eqn]$ -subproblem, it can be obtained that

$$[eqn]$$

It can be seen that

$$[eqn]$$

Considering the $[eqn]$ -subproblem, we obtain

$$[eqn]$$

and it implies that

$$[eqn]$$

Therefore, $[eqn]$ satisfies the following KKT conditions

$$[eqn]$$

The KKT conditions presented below can be applied to determine the termination criteria for Algorithm 1

$$[eqn]$$

where $[eqn]$ stands for a predefined tolerance. As has been mentioned previously, the sequence complies with the KKT conditions of the Lagrange function (25). $[eqn]$

Theorem 4. In our algorithm, the sequences $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ are Cauchy sequences and converge to their critical points.

Proof of Theorem 4: We first prove that $[eqn]$ is a Cauchy sequence. From the update rule $[eqn]$ , we can infer the following result

$$[eqn]$$

where $[eqn]$ .

Furthermore, this leads to the relation

$$[eqn]$$

where the first equality follows from $[eqn]$ , the second from Theorem 2, and the third from the unitary invariance of the Frobenius norm.

Recall the subproblem

$$[eqn]$$

The KKT conditions yield

$$[eqn]$$

where $[eqn]$ is the Lagrange multiplier.

If $[eqn]$ , it follow sthat

$$[eqn]$$

If $[eqn]$ and $[eqn]$ , it can be concluded that

$$[eqn]$$

If $[eqn]$ and $[eqn]$ , this leads to the conclusion that

$$[eqn]$$

Thus, in all cases, $[eqn]$ . Substituting into (94), we obtain the upper bound

$$[eqn]$$

Therefore,

$$[eqn]$$

Since $[eqn]$ is bounded by M1 and $[eqn]$ with $[eqn]$ and $[eqn]$ , it can therefore be concluded that

$$[eqn]$$

The right-hand side is a convergent geometric series, so $[eqn]$ is a Cauchy sequence.

Next, we prove that $[eqn]$ is a Cauchy sequence. Using the update $[eqn]$ , we have

$$[eqn]$$

Since $[eqn]$ is bounded and $[eqn]$ is a Cauchy sequence, it follows that

$$[eqn]$$

Hence, $[eqn]$ is a Cauchy sequence.

Subsequently, we provide a formal proof that $[eqn]$ is a Cauchy sequence. Similar to the proof of Lemma 9, we derive the following inequality.

$$[eqn]$$

Here, M2 is its upper bound. From the update $[eqn]$ , it implies

$$[eqn]$$

Then,

$$[eqn]$$

Utilizing the linearity property of the t-product, we can rewrite the expression as $[eqn]$ = $[eqn]$ . Within the framework of tensor t-product, the norm typically satisfies the following inequality

$$[eqn]$$

where $[eqn]$ denotes the spectral norm of tensor $[eqn]$ . Because $[eqn]$ is a fixed input tensor and its norm is bounded, then

$$[eqn]$$

Since it has been proven that $[eqn]$ is a Cauchy sequence, the right-hand side is bounded.

On the basis of (109) and $[eqn]$ is a bounded sequence, we have

$$[eqn]$$

Therefore, $[eqn]$ is a Cauchy sequence.

Next, we prove that $[eqn]$ is a Cauchy sequence. From the update $[eqn]$ , this yields

$$[eqn]$$

where $[eqn]$ .

The subproblem for $[eqn]$ is

$$[eqn]$$

By Lemma 5, the optimality condition is

$$[eqn]$$

Note that $[eqn]$ if and only if $[eqn]$ .

If $[eqn]$ , then

$$[eqn]$$

Since $[eqn]$ and $[eqn]$ , we have $[eqn]$ . Moreover, by Lemma 5, there exists a positive lower bound $[eqn]$ such that $[eqn]$ when $[eqn]$ , so

$$[eqn]$$

Let M3 be the upper bound in (115). Then for $[eqn]$ ,

$$[eqn]$$

If $[eqn]$ , then $[eqn]$ , and the inequality holds trivially.

Therefore,

$$[eqn]$$

Since $[eqn]$ is bounded by M4 and $[eqn]$ grows geometrically, we have

$$[eqn]$$

Thus, $[eqn]$ is a Cauchy sequence.

Finally, the proof that $[eqn]$ is a Cauchy sequence follows similarly from the update $[eqn]$ .

In conclusion, $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ are all Cauchy sequences and hence converge to the critical points of the objective function $[eqn]$ . $[eqn]$

Additionally, by running Algorithm 1 on the four datasets above, we practically validated the convergence of the proposed 2D-NLRSC. As shown in Fig 7, the error curves quickly stabilize after a few iterations on each dataset. These two aspects collectively demonstrate that the 2D-NLRSC model exhibits stability and rapidity.

6.1 Convergence rate analysis

We complement the convergence proof by quantifying the rates for both feasibility residuals and objective/value descent under standard assumptions used by augmented Lagrangian and splitting schemes.

Notation 1. Let the penalty parameters satisfy $[eqn]$ and $[eqn]$ with a fixed $[eqn]$ . Define the primal feasibility residuals

$$[eqn]$$

and the data-consistency residual

$$[eqn]$$

Theorem 5 (Geometric decay of feasibility residuals). Suppose $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ are bounded (as established by Lemmas 6–9). Then there exist finite constants $[eqn]$ such that, for all $[eqn]$ ,

$$[eqn]$$

Consequently, since $[eqn]$ and $[eqn]$ with $[eqn]$ , all residuals decay R-linearly

$$[eqn]$$

Proof: From the optimality of the $[eqn]$ -subproblem and the update $[eqn]$ , we have

$$[eqn]$$

The same argument using the $[eqn]$ and $[eqn]$ updates yields the bounds for $[eqn]$ and $[eqn]$ . For $[eqn]$ , the optimality of the $[eqn]$ -subproblem together with $[eqn]$ gives

$$[eqn]$$

Geometric decay follows from $[eqn]$ . $[eqn]$

Remark 1. Theorem 5 implies that the KKT termination test in (92) is met after $[eqn]$ iterations, dominated by the largest of the four residuals.

7 Conclusion

This study proposes 2D-NLRSC, a submodule clustering framework using nonconvex low-rank tensor approximation. Unlike traditional vectorization approaches, it rotates 2D images into tensors and employs the t-product to derive self-representative tensors. The key innovations involve an $[eqn]$ -induced tensor nuclear norm for optimal low-rank representation, along with a combined $[eqn]$ -norm and $[eqn]$ -norm regularization to enhance structure capture while effectively handling noise and redundancy. Theoretical convergence is demonstrated through the KKT conditions. Experiments on ORL, JAFFE, CMU-PIE, Yale, and MNIST show superior clustering accuracy and computational efficiency. Despite these advantages, this study still has certain limitations. The performance of 2D-NLRSC relies on the selection of several hyperparameters, which may require adjustment when applied to new datasets. Moreover, although relatively efficient, the computational complexity of its iterative optimization process remains relatively high for very large-scale image data. Future work will focus on automating hyperparameter selection and extending the framework to handle higher-order imaging data.