Reduced Basis Decomposition: a Certified and Fast Lossy Data Compression Algorithm

Dimension reduction is often needed in the area of data mining. The goal of these methods is to map the given high-dimensional data into a low-dimensional space preserving certain properties of the initial data. There are two kinds of techniques for this purpose. The first, projective methods, builds an explicit linear projection from the high-dimensional space to the low-dimensional one. On the other hand, the nonlinear methods utilizes nonlinear and implicit mapping between the two spaces. In both cases, the methods considered in literature have usually relied on computationally very intensive matrix factorizations, frequently the Singular Value Decomposition (SVD). The computational burden of SVD quickly renders these dimension reduction methods infeasible thanks to the ever-increasing sizes of the practical datasets. In this paper, we present a new decomposition strategy, Reduced Basis Decomposition (RBD), which is inspired by the Reduced Basis Method (RBM). Given $X$ the high-dimensional data, the method approximates it by $Y \, T (\approx X)$ with $Y$ being the low-dimensional surrogate and $T$ the transformation matrix. $Y$ is obtained through a greedy algorithm thus extremely efficient. In fact, it is significantly faster than SVD with comparable accuracy. $T$ can be computed on the fly. Moreover, unlike many compression algorithms, it easily finds the mapping for an arbitrary ``out-of-sample'' vector and it comes with an ``error indicator'' certifying the accuracy of the compression. Numerical results are shown validating these claims.


Introduction
Dimension reduction is ubiquitous in many areas ranging from pattern recognition, clustering, classification, to fast numerical simulation of complicated physical phenomena. The fundamental question to address is how to approximate a n-dimensional space by a d-dimensional one with d n. Specifically, we are given a set of high-dimensional data and the goal is to find its low-dimensional approximation Y = [y 1 , y 2 , . . . , y d ] ∈ R m×d (2) with reasonable accuracy. There are two types of dimension reduction methods. The first category consists of "projective" ones. These are the linear methods that are global in nature, and that explicitly transform the data matrix X into a low-dimensional one by Y = T X. The leading examples are the Principal Component Analysis (PCA) and its variants. The methods in the second category act locally and are inherently nonlinear. For each sample in the high-dimensional space (e.g. each column of X), they directly find their lowdimensional approximations by preserving certain locality or affinity between nearby points.
In this paper, inspired by the reduced basis method (RBM), we propose a linear method called "Reduced Basis Decomposition (RBD)". It is much faster than PCA/SVD-based techniques. Moreover, its low-dimensional vectors are equipped with error estimator indicating how close they are approximating the high-dimensional data. RBM is a relative recent approach to speed up the numerical simulation of parametric Partial Differential Equations (PDEs) [16,19,20,6,5]. It utilizes an Offline-Online computational decomposition strategy to produce surrogate solution (of dimension N ) in a time that is of orders of magnitude shorter than what is needed by the underlying numerical solver of dimension N N (called truth solver hereafter). The RBM relies on a projection onto a low dimensional space spanned by truth approximations at an optimally sampled set of parameter values [1,8,17,18,13]. This low-dimensional manifold is generated by a greedy algorithm making use of a rigorous a posteriori error bounds for the field variable and associated functional outputs of interest which also guarantees the fidelity of the surrogate solution in approximating the truth approximation.
The RBD method acts in a similar fashion. Given the data matrix X as in (1), it iteratively builds up Y (2) whose column space approximates that of X. It starts with a randomly selected column of X (or a user input if existent). At each step where we have k vectors {y 1 , . . . , y k }, the next vector y k+1 is found by scanning the columns of X and locating the one whose error of projection into the current space span{y 1 , . . . , y k } is the largest. This process is continued until the maximum projection/compression error is small enough or until the limit on the size of the reduced space is reached. An important feature is an offline-online decomposition that allows the computation of the compression error, and thus the cost of locating y k+1 , to be independent of (the potentially large) m.
This paper is organized as follows. In Section 2, we review the background material, mainly the RBM. Section 3 describes the reduced basis decomposition algorithm and discuss its properties. Numerical validations are presented in Section 4, and finally some concluding remarks are offered in Section 5.

Background
The reduced basis method was developed for use with finite element methods to numerically solve PDEs. We assume, for simplicity, that the problems (usually parametric partial differential equations (PDE)) to simulate are written in the weak form: find u(µ) in an Hilbert space X such that a(u(µ), v; µ) = f (v; µ), ∀v ∈ X where µ is an input parameter. These simulations need to be performed for many values of µ chosen in a given parameter set D. In this problem a and f are bilinear and linear forms, respectively, associated to the PDE (with a N and f N denoting their numerical counterparts). We assume that there is a numerical method to solve this problem and the solution u N , called the "truth approximation" or "snapshot", is accurate enough for all µ ∈ D.
The fundamental observation utilized by RBM is that the parameter dependent solution u N (µ) is not simply an arbitrary member of the infinitedimensional space associated with the PDE. Instead, the solution manifold M = {u N (µ), µ ∈ D} can typically be well approximated by a lowdimensional vector space. The idea is then to propose an approximation of M by W N = span{u N (µ 1 ), . . . , u N (µ N )} where, u N (µ 1 ), . . . , u N (µ N ) are N ( N ) pre-computed truth approximations corresponding to the parameters {µ 1 , . . . , µ N } judiciously selected according to a sampling strategy [13]. For a given µ, we now solve in W N for the reduced solution u (N ) (µ). The online computation is N -independent, thanks to the assumption that the (bi)linear forms are affine 2 and the fact that they can be approximated by affine (bi)linear forms when they are nonaffine [2,10]. Hence, the online part is very efficient. In order to be able to "optimally" find the N parameters and to assure the fidelity of the reduced basis solution u (N ) (µ) to approximate the truth solution u N (µ), we need an a posteriori error estima- and stability information of the bilinear form [12,14,19,20,21]. With this estimator, we can describe briefly the classical greedy algorithm used to find the N parameters µ 1 , . . . , µ N and the space W N . We first randomly select one parameter value and compute the associated truth approximation. Next, we scan the entire (discrete) parameter space and for each parameter in this space compute its RB approximation u (N =1) and the error estimator ∆ 1 (µ). The next parameter value we select, µ 2 , is the one corresponding to the largest error estimator. We then compute the truth approximation and thus have a new basis set consisting of two elements. This process is repeated until the maximum of the error estimators is sufficiently small.
The reduced basis method typically has exponential convergence with respect to the number of pre-computed solutions [15,4,3]. This means that the number of pre-computed solutions can be small, thus the computational cost reduced significantly, for the reduced basis solution to approximate the finite element solution reasonably well. The author and his collaborators showed [7] that it works well even for a complicated geometric electromagnetic scattering problem that efficiently reveals a very sensitive angle dependence (the object being stealthy with a particular configuration).

Reduced basis decomposition
In this section, we detail our proposed methodology by stating the algorithm, studying the error evaluation, and pinpointing the computational cost.

The algorithm
At the heart of the method stated in Algorithm 1 is a greedy algorithm similar to that used by RBM. It builds the reduced space dimension-bydimension. At each step, the greedy decision for the best next dimension to pursue in the space corresponding to the data is made by examining an error indicator quantifying the discrepancy between the uncompressed data and the one compressed into the current (reduced) space.
In the context of the RBM, we view each column (or row if we are compressing the row space) of the matrix as the fine solution of certain (virtual) parametric PDE with the (imaginary) parameter taking a particular value. Since this solution is explicitly given already by the data, the fact that the PDE and the parameter are absent does not matter. Once this common mechanism satisfied by each column (or row) is identified, the greedy algorithm still relies on an accurate and efficient estimate quantifying the error between the original data and the compressed one. This will be the topic of the next subsection.
To state the algorithm, we assume that we are given a data matrix X ∈ R m×n , the largest dimension d max < n that the practitioner wants to retain, and a tolerance R capping the discrepancy between the original and the compressed data. The output is the set of bases for the compressed data (a low-dimensional approximation of the original data) Y ∈ R m×d and the transformation matrix T ∈ R d×n . Here, d ≤ d max is the actual dimension of the compressed data.
With this output, we can Compress. We represent any data entry X(:, j), the j th column of X ∈ R m , by the j th column of T , T (:, j) ∈ R d , with usually d m.
Uncompress. An approximation of the data is reconstructed by , E cur = +∞, and i a random integer between 1 and n.

2.2.
Apply the modified Gram-Schmidt orthonormalization to obtain the d th basis of the compressed space Evaluate the compression of out-of-sample data. Given any v ∈ R m×1 that is not equal to any column of X, its compressed representation in

Efficient quantification of the error
A critical part to facilitate the greedy algorithm and make the algorithm realistic is an efficient mechanism measuring (or estimating) the error v − v C under certain norm, v − v C , in Step 2.3 of the algorithm. In this work, we are using the A−norm defined as follows. For a given symmetric and positive definite matrix A ∈ R m×m , the A−norm of a vector v ∈ R m×1 is defined by For v being any column of the data matrix X and v C its low-dimensional This makes the quality of compression uniform. In this case, the evaluation of (3) is greatly simplified and the algorithm is the fastest as shown below by the numerical results.
2. General diagonal matrix: This setting can be used if part of each data entry needs to be preserved better and other parts can afford less fidelity.
3. General SPD matrix: This most general case can be helpful if the goal is to preserve data across different entries anisotropiclly.
The goal is then to evaluate the error through (3) as efficiently as possible for any given c. This is achieved by employing an offline-online decomposition strategy where the c-independent parts are evaluated beforehand (offline) enabling a quick turnaround time for any given c encountered online. The specifics are given in the next subsection.

Computational cost and implementation aspects
The Offline-Online decomposition of the computations and their complexities are as follows. Here, we use nnz(A) to denote the number of nonzero entries of a sparse matrix A.
Offline The total cost is of order Offline MGS Every basis needs to orthogonalized against the current set of bases. The total cost is of order m d 2 max . Offline Calculation of Errors The next basis is located by comparing each column with its compressed version into the current space. To enable that, we encounter the following computational cost: Pre-computation of diag(X AX) (for v Av in (3)) and AX (for Av in (3)). The cost is of order nnz(A) n. Expansion of Y AX and Y AY . The former takes time of order nnz(A) d max m, and the latter of order nnz(A) d 2 max . Offline Searching After these calculations, the comparison between the original and compressed data is then only dependent on the size of c (which is also the number of columns for Y ). The complexity is of order d 2 max . It will be repeated for up to n times in the searching process of step 2.3 of the algorithm for each of the up to d max basis elements. The total cost is at the level of n d 3 max .
Online Given any (possibly out-of-sample) data v ∈ R m×1 , its coefficients in the compressed space is obtained by evaluating c = Y v. The cost is of order m d max . The decoding (Y c) can be done with the same cost.
The online computation has complexity of order m d max .
We remark that, if the actual practice does not requires forming v C (e.g. clustering and classification etc) and so we only work with the coordinates c of v in the compressed space, then the online cost will be independent of m and thus much smaller.

Numerical Results
In this section, we test the reduced basis decomposition on image compression, and data compression. Lastly, we devise a simple face recognition algorithm based on RBD and test it on a database of 575 images while comparing RBD with 6 other face recognition algorithms. The computation is done, and thus the speedup numbers reported herein should be understood as, in Matlab 2014a on a 2011 IMac with a 3.4 GHz Intel Core i7 processor.

Image Compression and comparison with SVD
We first test it on compressing two standard images Lena and Mandrill in Figure 1. They both have an original resolution of 512 × 512. We take A = I and test the algorithm. For each component of every image, we run the algorithm with d max ∈ {170, 51, 25} which implies a compression ratio of 33%, 10%, and 5% respectively. The resulting images (formed by multiplying the corresponding Y and T together) are shown on the 1 st and 3 rd row of Figures 2. As a comparison, we run SVD and obtain the reconstructed matrices with the first d max singular values accordingly. The resulting images are on the second and last row. Clearly, SVD provides the best quality pictures among all possible algorithms (and thus better than what RBD provides). However, we see that the RBD pictures are only slightly blurrier. Moreover, it takes much less time. In fact, we show the comparison in time between SVD and RBD in Table 1. We see that, when d = 51, RBD is three times faster than SVD and seven times faster when d = 25. Here the SVD time is the shorter between those taken by svd and svds commands in Matlab.

Data Compression
Here, we test the algorithm on a few artificially-generated data sets. Given a function f (x, y), the data denoted by f (D) is constructed by evaluating f on a uniform tensorial grid D := (x i , y j ) n i,j=1 .

Exact reconstruction
For tensorial functions such as those listed in Table 2 with their corresponding d values, the RBD method detects the optimal dimension, stops the greedy algorithm after d steps and decompose the matrix f (D) accordingly, that is, as an exact product of n × d and d × n matrices.

Face Recognition
Here, we demonstrate the superior efficiency and accuracy of the RBD method on a classical classification task -face recognition. The goal of face recognition is to recognize subjects based on facial images. It has important applications in areas ranging from surveillance, authentication, to humancomputer interaction etc. We use the UMIST database [9] that is publicly available on Roweis' web page 4 . Table 3 summarizes its characteristics: It contains 20 people under different poses. The number of different views per subject varies from 19 to 48. We use the cropped version whose snapshot is shown in Figure 4.
As in [11], we randomly choose 10 views from each class to form a training set. The rest of the samples (375 of them) are used as testing images. We show the average classification error rates in Figure 5 Left. These averages are computed over 100 random formations of the training and test sets. Shown in the middle are the results of six traditional dimension reduction techniques taken from [11]. Clearly, our method has similar performance as the PCA 9. Beyond spectral methods and trace optimization. While this paper focused on dimension reduction based on spectral techniques and trace optimization, other existing powerful methods rely on convex optimization with constraints. This section briefly describes two examples in this class for illustration purposes. For a recent survey of these techniques see [7] for example.
Possibly the best known technique along these lines in supervised learning is the method of Support Vector Machines (SVM); see [8,12,48].
It is in spirit similar to LDA (cf. Section 5.2) in that it finds a one dimensional projection to separate the data in some optimal way. Formally, the SVM approach consists of finding a hyperplane which best separates two training sets belonging to two classes. If the hyperplane is w T x + b = 0, then the classification function would be f (x) = sign(w T x + b). This will assign the value y = +1 to one class and y = −1 to the other, and it is capable of perfectly separating the two classes in ideal situations when the classes are linearly separable.
One of the key ingredients used by SVM is the notion of margin, which is the distance between two parallel support planes for the two classes. First, observe that the parameters w, b can be normalized by looking for hyperplanes of the form w T x + b ≥ 1 to include one set and w T x + b ≤ −1 to include the other. With y i = +1 for one class and y i = −1 for method, and outperforms three of the other five methods. However, RBD is much faster than PCA and other methods since they all involves solving eigenproblems [11]. A speedup factor as a function of the number of bases is plotted in Figure 5 Right which demonstrates a speedup factor of larger than two for this particular test when we reach the asymptotic region (around when the number of basis vectors is 25).

Concluding remarks
This paper presents and tests an extremely efficient dimension reduction algorithm for data processing. It is multiple times faster than the SVD/PCAbased algorithms. What makes this possible is a greedy algorithm that iteratively builds up the reduced space of basis vectors. Each time, the next dimension is located by exploring the errors of compression into the current space for all data entries. Thanks to an offline-online decomposition mechanism, this searching is independent of the size of each entry. Numerical results including one concerning a real world face recognition problem confirm these findings.