An Efficient Algorithm for Solving the Matrix Optimization Problem in the Unsupervised Feature Selection

Throughout this paper, we use $R^{n\times m}$ to denote the set of $m\times n$ real matrices. We write $B\ge 0$ if the matrix B is nonnegative. The symbols $Tr(B),B^T$ stand for the trace and transpose of the matrix B, respectively. The symbol $\parallel \alpha \parallel$ stands for the $l_2$-norm of the vector $\alpha$, i.e., $\parallel \alpha \parallel ={({\alpha }^T\alpha )}^{\frac{1}{2}}$. The symbol ${\parallel B\parallel }_F$ stands for the Frobenius norm of the matrix B. The symbol $I_p$ stands for the $p\times p$ identity matrix. For the matrices A and B, $A\odot B$ denotes the Hadamard product of A and B. The symbol $max\{x,y\}$ represents the greater of x and y.

In this paper, we consider the following symmetric matrix optimization problem in unsupervised feature selection.

Problem 1 arises in unsupervised feature selection, which is an important part of machine learning. This can be stated as follows. Data from image processing, pattern recognition and machine learning are usually high-dimensional data. If we deal with these data directly, this may increase the computational complexity and the memory of the algorithm. In particular, it may lead to the overfitting phenomenon for the machine learning model. Feature selection is a common dimension reduction method, the goal of which is to find the most representative feature subset from the original features, that is to say, for a given original high-dimensional data matrix A, we must find out the relationship between the original feature space and the subspace generated by the selected feature. Feature selection can be formalized as follows. where I denotes the index set of the selected features and Y is the coefficient matrix of the initial feature space in the selected features. From the viewpoint of matrix factorization, feature selection is expressed as follows Considering that the data in practical problems are often nonnegative, we add a constraint to guarantee that any feature is described as the positive linear combination of the selected features, so the problem in (3) can be rewritten as in (1).

In the last few years, many numerical methods have been proposed for solving optimization problems with nonnegative constraints, and these methods can be broadly classified into two categories: alternating gradient descent methods and alternating nonnegative least squares methods. The most commonly used alternating gradient descent method is the multiplicative update algorithm [1,2]. Although the multiplicative update algorithm is simple to implement, it lacks a convergence guarantee. The alternating nonnegative least squares method is used to solve nonnegative subproblems. Many numerical methods, such as the active method [3], the projected gradient method [4,5], the projected Barzilai–Borwein method [6,7], the projected Newton method [8] and the projected quasi-Newton method [9,10], have been designed to solve these subproblems.

For optimization problems with orthogonal constraints, which are also known as optimization problems on manifolds, there are many algorithms to solve this type of problem. In general, these can be divided into two categories: the feasible method and the infeasible method. The feasible method means that the variance obtained after each iteration must satisfy the orthogonal constraint. Many traditional optimization algorithms, such as the gradient method [11], the conjugate method [12], the trust-region method [13], Newton method [8] and the quasi-Newton method [14], can be used to deal with optimization problems on manifolds. Wen and Yin [15] proposed the CMBSS algorithm, which combined the Cayley transform and the curvilinear search approach with BB steps. However, the computational complexity increases when the number of variables or the amount of data increases, resulting in the low efficiency of this algorithm. Infeasible methods can overcome this disadvantage when facing high-dimensional data. In 2013, Lai and Osher [16] proposed a splitting method based on Bregman iteration and the ADMM method for orthogonality constraint problems. The SOC method is a valid and efficient method for solving the convex optimization problems but the proof of its convergence is still uncertain. Thus, Chen et al. [17] put forward a proximal alternating augmented Lagrangian method to solve such optimization problems with a non-smooth objective function and non-convex constraint.

Some unsupervised feature selection algorithms based on matrix decomposition have been proposed and have achieved good performance, such as SOCFS, MFFS, RUFSM, OPMF and so on. Based on the orthogonal basis clustering algorithm, SOCFS [18] does not explicitly use the pre-computed local structure information for data points represented as additional terms of their objective functions, but directly computes latent cluster information by means of the target matrix, conducting orthogonal basis clustering in a single unified term of the objective function. In 2017, Du S et al. [19] proposed RUFSM, in which robust discriminative feature selection and robust clustering are performed simultaneously under the $l_{2,1}$-norm, while the local manifold structures of the data are preserved. MFFS was developed from the viewpoint of subspace learning. That is, it treats feature selection as a matrix factorization problem and introduces an orthogonal constraint into its objective function to select the most informative features from high-dimensional data.

OPMF [20] incorporates matrix factorization, ordinal locality structure preservation and inner-product regularization into a unified framework, which can not only preserve the ordinal locality structure of the original data, but also achieve sparsity and low redundancy among features.

However, research into Problem 1 is very scarce as far as we know. The greatest difficulty is how to deal with the nonnegative and orthogonal constraints, because the problem has highly structured constraints. In this paper, we first use the penalty technique to deal with the orthogonal constraint and reformulate Problem 1 as a minimization problem with nonnegative constraints. Then, we design a new method for solving this problem. Based on the auxiliary function, the convergence theorem of the new method is derived. Finally, some numerical examples show that the new method is feasible. In particular, some simulation experiments in unsupervised feature selection illustrate that our algorithm is more efficient than the existing algorithms.

The rest of this paper is organized as follows. A new algorithm is proposed to solve Problem 1 in Section 2 and the convergence analysis is given in Section 3. In Section 4, some numerical examples are reported. Numerical tests on the proposed algorithm applied to unsupervised feature selection are also reported in that section.

In this section we first design a new algorithm for solving Problem 1; then, we present the properties of this algorithm.

Problem 1 is difficult to solve due to the orthogonal constraint $X^{T}X=I_p$. Fortunately, this difficulty can be overcome by adding a penalty term for the constraint. Therefore, Problem 1 can be transformed into the following form where $\rho$ is a penalty coefficient. This extra term is used to penalize the divergence between $X^{T}X$ and $I_p$. The parameter $\rho$ is chosen by users to make a trade-off between making $\frac{1}{2}{\Vert A-AXY\Vert }_F^2$ small, while ensuring that $\Vert X^{T}X-I_p{\Vert }_F^2$ is not excessively large.

Let the Lagrange function of (4) be where $\alpha \in R^{m\times p}$ and $\beta \in R^{p\times m}$ are the Lagrangian multipliers of X and Y. It is straightforward to obtain its gradient functions as follows Setting the partial derivatives of X and Y to zero, we obtain which implies that Noting that $(X,Y)$ is the stationary point of (4), if it satisfies the KKT conditions which implies or According to (6) and (7) we can obtain the following iterations which are equivalent to the following update formulations However, the iterative formulae (10) and (11) have two drawbacks, as follows.

In order to overcome these difficulties, we designed the following iterative methods where Here, $\sigma$ is a small positive number which can guarantee the nonnegativity of every element of X and Y. So we can establish a new algorithm for solving Problem 1 as follows (Algorithm 1).

Algorithm 1: This Algorithm attempts to Solve Problem 1.

Input Data matrix $A\in R^{n\times m}$, the number of selected features p, parameters $\sigma$, $\rho$ and $\delta$. Output An index set of selected features $I\subseteq \{1,2,\cdots ,m\}$ and $\left|I\right|=p.$ 1. Initialize matrix $X^{(0)}\ge 0$ and $Y^{(0)}\ge 0$. 2. Set $k:=0$. 3. Repeat 4. Fix Y and update X by (2.9); 5. Fix X and update Y by (2.10); 6. Until convergence condition has been satisfied, otherwise set $k:=k+1$ and turn to step 3. 7. End for 8. $X={(x_1,x_2,\dots ,x_n)}^T.$ Compute $\Vert x_i\Vert$ and sort them in a descending order to choose    the top p features.

The sequences $\left\{X^{(k)}\right\}$ and $\left\{Y^{(k)}\right\}$ generated by Algorithm 1 have the following property.

Case I.

From ${[{\nabla }_{X}F(X,Y)]}_{ij}\ge 0$ it follows that ${\overline{X}}_{ij}=X_{ij}$ and

Hence, if $X_{ij}^{(k)}>0$ then $X_{ij}^{(k+1)}>0$ and if $X_{ij}^{(k)}\ge 0$ then $X_{ij}^{(k+1)}\ge 0$.

Case II.

From ${[{\nabla }_{X}F(X,Y)]}_{ij}<0$, it follows that ${\overline{X}}_{ij}\ne X_{ij}$ and

Noting that $max\{X_{ij}^{(k)},\sigma \}>0$ and ${[{\nabla }_{X}F(X^{(k)},Y^{(k)})]}_{ij}<0$, we can easily conclude that if $X_{ij}^{(k)}>0$ then $X_{ij}^{(k+1)}>0$ and if $X_{ij}^{(k)}\ge 0$ then $X_{ij}^{(k+1)}\ge 0$.

Case III.

From ${[{\nabla }_{Y}F(X,Y)]}_{ij}\ge 0$ it follows that ${\overline{Y}}_{ij}=Y_{ij}$ and

Thus, if $Y_{ij}^{(k)}>0$ then $Y_{ij}^{(k+1)}>0$ and if $Y_{ij}^{(k)}\ge 0$ then $Y_{ij}^{(k+1)}\ge 0$.

Case IV.

From ${[{\nabla }_{Y}F(X,Y)]}_{ij}<0$, it follows that ${\overline{Y}}_{ij}\ne Y_{ij}$ and

Based on the fact that $max\{Y_{ij}^{(k)},\sigma \}>0$ and ${[{\nabla }_{Y}F(X,Y)]}_{ij}<0$, we can conclude that if $Y_{ij}^{(k)}>0$ then $Y_{ij}^{(k+1)}>0$ and if $Y_{ij}^{(k)}\ge 0$ then $Y_{ij}^{(k+1)}\ge 0$.  □

In this section, we will give the convergence theorem for Algorithm 1. For the objective function of Problem (4), we first prove that where $X^{(k)}$ and $Y^{(k)}$ are the k-th iteration of Algorithm 1, then obtain the limit point as the stationary point of Problem (4). In order to develop this section, we need a lemma.

Similarly, we can use the same method to verify that when Y is fixed, the function $F(X,Y)$ is also a non-increasing function. Thus, we have Consequently, by (20) and (21), we obtain

In this section, we first present a simple example to illustrate that Algorithm 1 is feasible to solve Problem 1, and we apply Algorithm 1 to unsupervised feature selection. We also compare our algorithm with the MaxVar Algorithm [22], the UDFS Algorithm [23] and the MFFS Algorithm [24]. All experiments were performed in MATLAB R2014a on a PC with an Intel Core i5 processor at 2.50 GHz with a precision of $\epsilon =2.22\times {10}^{-16}$. Set the gradient value Due to the KKT conditions (5) we know that if $GV(X,Y)=0$ then $(X,Y)$ is the stationary point of Problem 1. So we use either $GV(X,Y)\le 1.0\times {10}^{-4}$ or the iteration step k has reached the upper limit 500 as the stopping criterion of Algorithm 1.