Hybrid Matrix Completion Model for improved Images Recovery and Recommendation Systems

Matrix completion methods have been widely applied in images recovery and recommendation systems. Most of them are only based on the low-rank characteristics of matrices to predict the missing entries. However, these methods lack consideration of local information. To further improve the performance of matrix completion. In this paper, we propose a novel model based on matrix decompositions and matrix local information. Specifically, we update a number of rank-one matrices, which circumvented the rank estimation in matrix decomposition. And a penalty function is designed to punish singular values without introducing additional parameters. The local information component extracts similar information by an adaptive filter via convolution operation which kernel is obtained by the minimum variance. Finally, we integrate matrix decomposition and local information components via different weights. We apply the proposed method to real-world image datasets and recommendation system datasets. The experimental results demonstrate the proposed model has a lower error and better robustness than several competing matrix completion methods.


I. INTRODUCTION
Matrix completion is to restore the missing entries in a sparse matrix based on partial observation entries, and has been enjoyed widespread applications in recent years, such as recommendation systems [1], [2], image processing [3]- [5], social networks [6], [7], large scale classification [8] and clustering [9]. In these applications, authors obtain the completed matrix by assuming the observed matrix has a lowrank or approximately low-rank structure. Thus, they use these methods, such as rank minimization (RM), nuclear norm minimization (NNM) and matrix factorization (MF), to predict the missing entries. The NNM is a relaxation to RM, and the MF decomposes a matrix into several factor matrices. Let Ω denotes the set of observed indices in an observation matrix, is a projection operator onto Ω , the rank minimization method is where ∈ ℝ × is a sparse observed matrix, is a completion matrix, i.e. deserved matrix.
( ) is the rank of . Model (1) completes the observed matrix by minimizing ( ). However, model (1) was often circumvented as an NNM due to the nonconvexity of rank function [10]. Candès et al. [11] proved that the RM problem can be surrogated by NNM. Thus, model (1) can be rewritten as follows, min ‖ ‖ * s. t.
where ‖ ‖ * is the nuclear norm of . At present, NNM has achieved satisfactory performance on matrix completion problems. For example, Cai et al [12] proposed a singular value thresholding algorithm (SVT) that using a softthresholding rule to punish singular values, constrain the nuclear norm and limit the rank size. However, the SVT suffers from high computational costs because of singular value decomposition (SVD) in each iteration. To enhance the efficiency of SVD in matrix completion, Oh et al. [13] only compute the partial singular values from a small factored matrix to avoid computation of SVD, which can improve the speed of SVT. The nuclear norm is the ℓ -norm of the singular vector according to the definition. And the nuclear norm treats every singular value as equally important. Therefore, some studies have proposed the weighted nuclear norm minimization (WNNM) [14], [15], which weights the singular values and provides a large advantage over the nuclear norm. But WNNM introduces some additional parameters in the regularizer. Nie et al. [16] proposed a new matrix completion method based on non-convex relaxation, using a non-convex function to punish the singular values without additional parameters. Besides, the NNM problems can be surrogated by the truncated nuclear norm minimization [17] and scaled latent nuclear norm minimization [18]. However, these NNM-based methods lack consideration of the local information when matrix completion. MF is another matrix completion method that treats the deserved matrix as a product of multiple factor matrices. For example, the deserved matrix is decomposed into ∈ ℝ × and ∈ ℝ × , satisfying = , where r is the rank of . The optimal rank will be determined by adjusting r. Besides, can also be decomposed into ∈ ℝ × , ∈ ℝ × and ∈ ℝ × . That is, = .
Where is a diagonal matrix, the values on its diagonal are singular values, the size of r can be restricted by changing the number of elements on the diagonal of . Since the sum of the values on the diagonal of is the nuclear norm, the WNNM can also be weighted to the diagonal of . The MF-based methods explore the top r singular values and require fewer parameters tuning, which is simpler to implement and use.
In the MF, the low-rank matrix fitting algorithm (LMaFit) is one classic matrix decomposition model that estimates the rank by two factor matrices [19]. Inspired by LMaFit method, Tanner et al. [20] proposed the alternating steepest descent algorithm (ASD) which can replace the least square subproblem solution in LMaFit. However, the ASD is a linear method for matrix completion. In order to handle nonlinear structure data, Fan et al. [21] proposed a deep matrix factorization (DMF) which uses a deep-structure neural network to complete a sparse observation matrix. Matrix completion in a non-linear way also includes Riemannian optimization [5], [22]. SVD is often involved in matrix factorization, and the computation of SVD is burdensome work [24]. Yao et al. [25] proposed an accelerated and inexact soft-impute algorithm (AIS) for large-scale matrix and tensor completion. They use the power method [26] to accelerate calculation, and has faster convergence without SVD.
Currently, side information is used in matrix completion because it exploits the internal information of the observation matrix. Xu et al. [27] assume the column and row vectors in sparse matrix lie in the subspace spanned by the column and row vectors in ∈ ℝ × and ∈ ℝ × , respectively. and are side information matrices. The side information is introduced into the matrix completion and the completed problem is transformed into solving min ‖ ‖ * , where ∈ ℝ × , and ( ) = = ( ) . Chiang et al. [28] consider the robustness of matrix completion model, they assumed that the most given side information has noisy, so the residual matrix is introduced to satisfy = + , where can capture the information that side information fails to describe. Besides, other interests integrate the prior knowledge or data attributes to the factorization model to improve the performance of matrix completion [29], [30]. Different from side information, local information learns local characteristics from the perspective of local parts, which extracting information has more explanatory and more specific.
Although the aforementioned methods are successful in many areas, most of them simply ignore local information. Therefore, we propose a novel matrix completion method which based on the low-rank matrix factorization and local information filtering, named LAMC. First, we construct a lowrank matrix factorization model, which completes the observation matrix by iterating the rank-one matrix. The power method is used to speed up and a non-convex function is used to penalize the maximum singular value in the iterative process. Second, we construct an adaptive filter to complete missing entries. The optimal filter kernels are obtained by the minimum variances of the missing entries area. A convolution operation is rented to obtain approximations of missing entries. Finally, we integrate the results of low-rank matrix factorization and the convolution operation to complete the missing entries. Our contributions are summarized as follows: 1) The proposed model combining low-rank matrix decomposition and adaptive local filtering, which can capture more features of the observation matrix. Because the low-rank matrix decomposition and adaptive local filtering are based on the row-column relationship and local information to complete the missing entries, respectively. 2) The low-rank matrix factorization component updates the deserved matrix by rank-one matrices. And a non-convex function that punishes the singular values in each iteration. Different from WNNM, the non-convex function was designed does not introduce additional parameters.
3) The convolution operations are executed with different kernel matrices, which cleverly solved the problem of hyperparameter learning (best kernel size), and fully considered the correlation of local information in observation matrices. 4) A heuristic algorithm is used to better blend our components. And experiments verify that the proposed method has a better performance when applied to image recovery and recommendation systems. The outline of this paper is organized as follows. Section 2 describes the related work; Section 3 first presents the lowrank factorization component and then introduces the convolution operation using adaptive filters, after that, the ensemble module is presented. Experimental results on realworld image datasets and recommendation system datasets are performed in section 4. Section 5 ends with some conclusions and future research lines.

II. RELATED WORK
In this section we will introduce one of the matrix factorization methods, rank-one matrix completion (R1MC), as well as some common filters, which are related to our next work.

A. RANK-ONE MATRIX COMPLETION
R1MC is one of matrix factorizations which estimates a matrix into a sum of rank-one or rank-R matrices. The sparse observation matrix can be completed by the following model: where is the completed matrix, is the sum of R rank-one matrices. R is the rank with the best matrix completion result. The missing entries are updated by ( ), Ω is the set outside Ω.
and are column vectors. We summarized Algorithm 1 to solve Model (3).
In Algorithm 1, consists of R , subject to = 1,2, … , . is a diagonal matrix, = ( , , … , ) and consists of R . ‖ ‖ is the Frobenius norm of matrix . The advantages of the R1MC algorithm are that only one parameter (i.e. rank R) needs to be adjusted, and no step-size parameter setting requirement [31], [32]. However, (i). it needs to perform SVD (see step 4 in Algorithm 1), and the SVD is a time-consuming job. (ii). Each update needs to use the entries of in Ω . If the obtained entries of ( ) have a larger error in the previous iteration, it will lead to an increased error in the next iteration.

B. COMMON FILTERS
Filters are usually used in image processing for noise suppression, smoothing and filtering [33]. Commonly filters include: mean filter, Gaussian filter and median filter. The mean filter performs average smoothing on an image (matrix or tensor), with the advantages of high efficiency and simplicity. However, calculating the average value will blur the edge information and loss of important local features in an image.
The Gaussian filter uses the distribution of two-dimensional Gaussian function to smooth images, which replaces the pixel value of the point with the weighted average of the pixel neighborhoods, and the weight of each neighborhood pixel is dependent on the distance from the center point. The twodimensional Gaussian function is rotationally symmetric and has the same degree of smoothness in all directions. It does not change the edge direction of the original image, and guarantees the characteristics of feature points and edges.
The basic idea of the median filter is to replace the missing values with the median of the neighborhood values. This method can preserve edge details while removing noise. The median filter overcomes the problem of blurring of image details caused by the mean filter to a certain extent.
Filters have the advantage of extracting local information, but traditional filters set a fixed-size filter for an image (matrix or tensor), assuming that the relative importance of local information is consistent. This is unreasonable because in the real world, the size of areas with locally relevant information is not fixed. In addition, the filter matrix size is a hyperparameter, and adjustment of hyperparameters is not easy in machine learning. However, the adaptive filters can handle this problem well.
Matrix completion can be regarded as a signal processing of observation entries to obtain the missing entries. Therefore, filtering the observation entries can complete a sparse matrix. Unlike filters that fill in missing entries from the perspective of local information, low-rank matrix completion predicts missing entries from the low-rank characteristics. Therefore, combining two components can help to achieve a better result.

A. MODEL ARCHITECTURE
The proposed model is based on two perspectives. (i). Some sparse matrices satisfy the low-rank characteristics, and the missing entries can be predicted by the matrix factorization. (ii). There is a certain relationship between the local entries in the sparse matrix, and the missing entries can be completed by the convolution operation. Therefore, we construct an ensemble model to predict the missing entries.
The proposed model is summarized in Fig.1. The completion objects are incomplete images and recommendation systems datasets. The sparse observation matrix can be obtained from these datasets. Then a low-rank matrix factorization component updates the matrix to obtain the one recovery result via min ‖ ( ) − ( )‖ , and is a sum of multiple rank-one matrices. The adaptive filtering component obtains the best filter for each missing entry, the number of the missing entries is |Ω |, so we need |Ω | filters. After that, the convolution operation is performed on the best filters and associated missing entry areas to obtain . The final completion matrix will be obtained by weighting and . In application, we can directly use to display the recovered images or sort the rows of to obtain the recommendation list. The details are shown in the subsequent sections.

B. RANK-ONE MATRIX COMPLETION
According to the SVD, any matrix can be represented as a sum of a number of rank-one matrices and its singular vectors are orthogonal. However, a matrix can also be the sum of more rank-one matrices if the singular vectors are not enforced to be orthogonal. So we construct a new matrix completion model based on R1MC as follows, rating11 rating12 ··· rating1n rating21 rating22 ··· rating2n · · · · · · · · · · · · ratingm1 ratingm2 ··· ratingmn Xm2 ··· Xmn · · · · · · · · · X rank11 rank12 ··· rank1n rank21 rank22 ··· rank2n · · · · · · · · · · · · rankm1 rankm2 ··· rankmn Model (4) complete the missing entries in by ( ), and are the left and right rank-one vectors. And the left vectors and are not enforced to be orthogonal, the right vectors are the same. Thus, N is the number of rank-one matrices rather than the optimal rank. The (·) is a penalty function for singular values. Model (4) wants to use N rankone matrices to obtain matrix with minimization ‖ ( ) − ( )‖ . We use Algorithm 2 to solve Model (4), the following is the summary. In step 4 of Algorithm 2, The power method is used to calculate , and , because the largest singular value and corresponding leading left and right singular vectors can be computed efficiently by the power method [26]. We chose of largest singular value to update because consists of (∑ ( ) ), and we can optimize it step by step. The function (·) in step 5 of algorithm 2 is a penalty function that is vital for solving mode (4) because it governs the lowrank constraint or suppresses singular values. The traditional methods are to treat (·) as a soft or hard threshold operation [34]. These approaches extremely suppress small singular values because they penalize smaller singular values to zero. However, considering that large singular values usually represent the signal itself, small singular values are likely to be composed of noise, therefore, (i) large singular values should be treated similarly, while (ii) small singular values should be suppressed excessively. In addition, 0 ≤ ( ) ≤ . Therefore, we designed the following parameterless penalty function to suppress singular values:

C. PROOF ITERATIVE CONVERGENCE OF LOW RANK MATRIX FICTORIZATION
To verify the convergence of Algorithm 2 when solving model (4), the convergence analysis is given below. We transform the proof of iterative convergence into the following question: Proof: For algorithm 2, we let the SVD of ( ) − ( ) = ∑ , where r is the rank of ( ) − ( ), ≥ ≥ ⋯ ≥ > 0. Then, In step 4 of algorithm 2, we select the maximum singular value and leading left and right singular vectors to update = + . Thus, we obtain According to the step 5 and 6 of algorithm 2, we have = ( ), = and = . The projection operation is equal to 0, thus ( − ) = − . Then formula (8) can be rewritten into the following formula: Obviously

D. ADAPTIVE LOCAL FILTERING MODULE
In this section, we introduce the convolution operation to extract local features and how to construct an adaptive filter. For the sparse observation matrix , the kernel matrix is ∈ ( where is a positive integer and 2 + 1 ≤ min ( + , + ) , where padding is the extended size that makes one missing entry at the center of the filter, The value of padding is usually taken as missing, i.e. 0. The missing entry , = 0 can be completed to obtain , . The mean and Gaussian filter formula for obtain , is = ( ⊙ ′( − : + , − : + )). (10) Where is the best filter for the missing entry .
(·) is a function that sums all the entries in a matrix. ′( : , : ) is to select the area where the abscissa is a to b, and the ordinate is c to d. ′ is the padding matrix of . The value in the kernel matrix is set to in the mean filter, where (•) is the indicator function. For Gaussian filter, is normalized according to the Gaussian distribution. However, for the median filter, the values of are all set to 1, and then the Hadamard product is performed with ( − : + , − : + ) and the (·) is a function that selects the median in a sparse matrix, the formula is as follows. In the next, we will obtain the best filter kernel size for all missing entries. Traditional filtering methods often set a fixed-size kernel matrix [35], [36]. However, the fixed-size kernel has some shortcomings. (i) It does not take into account local information consistency. Because the fixed-size kernel sees all areas equally. (ii) It requires learning hyperparameters. The kernel size of the filter (2k+1) is a hyperparameter that requires expert experience to set. This is a tedious job and does not guarantee absolute reliability.
For formulas (10) and (11), it can be known that the essence of filtering is to obtain new values based on the local correlations in an observation matrix. To determine the relevance of a specific region in it, we calculate the variance of different ranges in each region, and the most relevant area corresponds to the smallest variance. Then we obtain where is the mean of ′( − : + , − : + ),. The value of formula (12) can reflect the correlation degree of region in ′( − : + , − : + ) . Fig.2(a) randomly removes 50% of the pixels and calculates the variance corresponding to different (see Fig.2(b)) while ( , ) = 0. The smallest can be determined when = 1, 2, … , min ( , ). That is, since has a certain dimension, we can definitely find a minimum variance corresponding to .

E. Ensemble learning module
In the above part, we obtain two completed matrices and by low-rank matrix factorization and adaptive local filtering. We hope that the completed matrix can be optimal as a whole, so we circumvent this ensemble problem by parameter learning and define the optimal weights and for and . The optimal weights are tuned through the differential evolution (DE) algorithm [37] because the DE algorithm has the advantage of quickly approaching the global optimal solution. DE algorithm includes three basic steps: mutation, crossover and selection. The individual differences in the population are large when the algorithm starts to iterate, and this mutation operation will make the algorithm have strong global search ability; the individual differences in the population are small, the algorithm to have a strong local search ability. The main advantages of DE are: (i). fewer parameters to be tuned, (ii). Not easy to fall into the local optimum, (iii). And faster convergence speed. The optimal weight vector is = ( , ) . Algorithm 3 shows the Variance process of selecting the optimal weights. The fitness function is defined as the smallest root square error.
Algorithm 3 Select the optimal weight through the DM 1: input: matrix X, Y; 2: initialization: population size NP, mutation factor F, crossover probability CR, maximum number of functions evaluations G, current number of iteration g = 0; 3: GENERATE random weight matrix , = 1,2, … , ; 4: while g < G do 5: MUTATION, generate mutation weight matrix = + ( − ), ≠ 1 ≠ 2 ≠ 3, 1, 2, 3 ∈ [1, ]; 6: CROSSOVER, The mutation weight matrix and the original weight matrix are crossed to produce the trial weight: , , ℎ , where j is the index of the j th weight, j = 1,2; 7: Selection, choose and with the smallest fitness value as the new weight , and g=g+1; 8: output: is the weight corresponding to the minimum fitness function value of ;
where Ω is the set of testing data set indices in sparse matrices. |Ω | denotes the number of entries in a testing set. All datasets are standardized so that the results of the low-rank matrix factorization component and the adaptive filtering component are in the same dimension.

A. GRAYSCALE IMAGES RECOVERY
In this part, we obtain three grayscale images from references [16] and [32], to compare the performance of our algorithm with other algorithms. The original images (see Fig.3) are named Lenna (512×512 pixels), Rice (600×960 pixels), and Male (1,024×1,024 pixels). 50% of each image data matrix entries are randomly sampled as the training set, another 5% as the validation set for tuning parameters, and the rest as the resting set. In our method, we apply the traditional filters to our adaptive filtering component. Table I shows the results of various filters applied to grayscale images. We can conclude from Table I: The adaptive mean filter has the best recovery performance for grayscale images. The adaptive median filter is the second most effective. The adaptive Gaussian filter has the worst effect on image recovery. Therefore, a combination of adaptive mean filter and lowrank matrix completion is selected to restore the missing entries in images. Table II shows the comparison between our proposed method and other matrix completion methods. According to table II, the LAMC algorithm performs best on each image data set (bold in the table). The overall performance of the AIS is second only to our method. The recovery performance of LMaFit and SVT are slightly similar, but the SVT is better than the LMaFit in general, for example, SVT performs better on the Rice and Lenna datasets. The performance of LDMM is unsatisfactory, especially on the Rice data set, which results are poor due to the limitation of the validation set. ASD has a poor performance on the three image datasets, and the error on the Male data set is larger than other algorithms. This group of experiments shows that our method has the best effect on grayscale image recovery. In order to verify the performance of different scale validation sets on the experiments, we keep the proportion of the training sets unchanged, set the proportion of the verification sets to 25%, and the test sets to 25%. Table III shows the recovery of the images through adaptive mean filtering and low-rank matrix completion. From the experimental results, it can be seen that the RMSE and MAE performance of LAMC are the best. The AIS algorithm has the second best performance in terms of image recovery, with errors that are only higher overall than those of our algorithm.

Lenna Rice Male
The SVT can also achieve better results on grayscale image recovery, although it is inferior to the ALS algorithm. The ASD has poor results on this proportion of datasets. The results obtained by the LDMM on the Rice image are more reliable due to the proportion increase of the validation sets. This experiment shows that our proposed algorithm has strong robustness and better performance. To fully demonstrate the completion results on images, we choose image datasets with a verification set ratio of 25% to show the recovery effect. Fig.4 shows the original sparse images, and Fig.5 shows the recovery images with different algorithms. It can be seen from Fig.5 that the LAMC has the best image restoration performance visually. The restored images obtained by the ASD are almost indistinguishable, and its restoration performance is the worst. The AIS and SVT can restore images better, and the restoration performance is only inferior to our proposed method from the visual comparison. The visual display of the image restoration matches with the RMSE of Table III.

Lenna
Rice Male  The low-rank matrix completion has been widely used in recommendation systems. However, the recommendation systems based on the matrix completion only consider the lowrank characteristic (i.e. the rows or columns of sparse observation matrix), in order to further improve the recommendation efficiency. We use common datasets, MovieLens 100K, MovieLens 1M [38] and Book-Crossing [39]. The MovieLens 100K contains 100,000 ratings issued by 943 users for 1,682 movies. The MovieLens 1M dataset consists of 1,000,209 ratings from 6,040 users on 3,703 movies. The Book-Crossing contains 13,123 users and 7,774 books, for a total of 1,048,574 user-book ratings.
We preprocess the datasets of the recommendation systems before conducting experiments. Because the random construction of the user-item rating matrices did not consider the local relations. Therefore, we defined a processing rule to put the most similar rows and the most similar columns together, that is, the upper row and the next row (the upper column and the next column) of each row (column) are the most similar rows (columns). The matrix constructed in this way helps to capture some local information and does not affect the row-column relations. All experiments select 60% of the rating data as the training sets, the proportion of the validation sets include 5% and 25% respectively, and the remaining parts are the test sets. Table IV shows the experimental results of different adaptive filter methods on the recommendation systems. It can be seen that the adaptive Gaussian filter performs best on the three recommendation system datasets, especially on the Book-Crossing data set. Furthermore, the adaptive Gaussian filter has the best RMSE on MovieLens 100K and Book-Crossing. However, the adaptive median filter has the best MAE on MovieLens 1M. The adaptive mean filter performs slightly inferior to the Gaussian filter on the MovieLens datasets, but the performance on the Book-Crossing dataset is not satisfactory. Experiments show that the adaptive Gaussian filter is more efficient for MovieLens datasets.  Table V shows the experimental results of adaptive Gaussian filter and low-rank matrix completion on the recommendation system with 5% validation sets. The RMSE and MAE of the proposed LAMC algorithm are smaller than other methods in general, especially for MovieLens 1M and Book-Crossing datasets. The SVT algorithm performs well on the MovieLens 100K, but it did perform well on MovieLens 1M. AIS and DMF perform better on the recommendation systems, and DMF has the best MAE on the MovieLens 1M. The LDMM cannot be performed on the Book-Crossing data set due to computer memory. ASD does not converge well on the Book-Crossing.
We integrate the adaptive Gaussian filter and low-rank matrix factorization to 25% of validation sets. Table VI displays the experimental results of the hybrid. The LAMC has the best RMSE on all recommendation system datasets, and its MAE is better in general. The performance of the DMF is second only to our method, and the MAE of DMF is the best on MovieLens 1M. The performance of the AIS is relatively stable, and will not cause large deviations with different proportions of validation sets. It demonstrates that AIS has better adaptability just like our method. The LDMM cannot be run on the Book-Crossing dataset due to computer memory. ASD does not converge well on the Book-Crossing data set. Experiments show the superiority of our method. To further verify the completion preference of our method, we show the RMSE versus time obtained from the low-rank matrix factorization. The MovieLens and Book-Crossing datasets with validation sets ratio of 25% are selected for this experiment. The experimental results are shown in Fig.6. It can be seen from fig.6 that the LAMC has better convergence and lower RMSE. The AIS has the fastest convergence on MovieLens 100K and Book-crossing datasets, but its RMSE is not better than ours. The DMF has the same performance as our method when we only perform matrix low-rank factorization, especially in the Book-Crossing data set, and its convergence speed is similar to ours. Because the LDMM cannot run on our machine on the Book-Crossing data set. Therefore, Fig.6(c) does not show the performance of the LDMM. The RMSE of the LMaFit fluctuates on the MovieLens 1M and Book-Crossing datasets, and fluctuates more obviously on the MovieLens 1M dataset. The ASD cannot achieve a good result within the specified time under the Book-Crossing data set, and its RMSE is more than 2.0. These experimental results show that our proposed method keeps a better convergence and lower RMSE without adaptive filtering. VOLUME XX, 2021

V. CONCLUSIONS
In this article, we propose a new matrix completion method that is a hybrid of low-rank matrix factorization and adaptive filtering. In the component of low-rank matrix factorization, the deserved matrix is decomposed as several rank-one matrices and then solved iteratively. The non-convex function and power method are used to punish singular values, reduce parameters and speed up computation. The convergence of the proposed low-rank matrix factorization is strictly proved mathematically. The second component of adaptive filtering which borrowing the idea of convolution to extract local information. We design a method to obtain the best kernel size for each missing entry. A heuristic algorithm is used to integrate the two components. We use real-world image datasets and recommendation system datasets to perform the experiments. And the results demonstrate that the proposed method is more accurate, much faster and more robust than the competing methods.
Since tensor is a high-dimensional matrix, we are interested in using an ensemble idea to high-dimensional tensor in future work. The tensor has a regular internal structure, and how to complete it without destroying the structure is also a matter we need to study. Besides, exploring the performance of different tensor norms to apply to high-dimensional, for example, color images and videos, etc., is one of our future works.