Interpretable Matrix Completion: A Discrete Optimization Approach

2018·arXiv

Abstract

1. Introduction

Low-rank matrix completion has attracted much attention after the successful application in the Netflix Competition. It is now widely utilized in far-reaching areas such as computer vision (Candes and Plan (2010)), signal processing (Ji et al. (2010)), and control theory (Boyd et al. (1994)) to generate a completed matrix from partially observed entries.

2 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

The classical low-rank matrix completion problem considers the following problem: Given a matrix with entries only partially known (denote Ω as the set of known entries), we aim to recover a matrix of rank k that minimizes a certain distance metric between X and A on the known entries of A:

where we normalized the objective so that it is O(1). The rank k constraint on X can be equivalently formulated as the existence of two matrices such that . Therefore, the problem can be restated as:

In many applications for matrix completion, it is customary for each row of the data to represent an individual and each column a product or an item of interest, and being the response data of individual i on item j. Therefore, the matrices U and V are commonly interpreted as the “user matrix” and “product matrix” respectively.

Let us denote each row of U as and each column as (similarly for V ). Then, () represents a “latent feature” for users (products), and in total there are k latent features for users (products). The goal of matrix completion is thus to discover such latent features of the users and the products, so that the dot product of such features on user i and product , is the intended response of user i on product j.

While this interpretation is intuitive, it does not offer insight on what the latent features of users and products mean. Inductive Matrix Completion, first considered in Dhillon et al. (2013), aims to rectify such problem by asserting that each of the k latent features is a linear combination of p > k known features. We focus on the one-sided information case, where only one of the user/product matrix is subject to such constraint. This case is more relevant as features related to users are often fragmented and increasingly constrained by data privacy regulations, while features about

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 3

products are easy to obtain. We would also have a short discussion later on why the two-sided information case is not interesting under the context of this paper.

Without loss of generality, we would assume the information is on the product matrix and denote the known feature matrix with p features as . As an example, if the items are movies, then represents feature j for a movie (actors, budget, running time, etc), and each product feature needs to be linear combination of such features. Mathematically, this translates to the constraint:

where . Therefore, the inductive version of the problem in (1) can be written as:

Although the inductive version of the problem adds more interpretability to the product latent features, they are still far from fully interpretable. For example, if we take the items to be movies, and features to be running time, budget, box office, and number of top 100 actors, then the generated features could look like:

These features, although a linear combination of interpretable features, are not very interpretable itself due to the involvement of multiple factors with different units. Therefore, it cannot signifi-cantly help decision makers to understand “what is important” about the product. Furthermore, the appearance of the same factor in multiple features with different signs (as shown above) further complicates any attempt at understanding the result.

Therefore, we argue that instead of supposing the k features are linear combinations of the p known features, we should assume that the k features are selected from the p known features. This

4 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

formulation alleviates the two previous problems mentioned: it guarantees the latent features to be interpretable (as long as the original features are), and it prevents any duplicating features in the selected k latent features. We denote this Interpretable Matrix Completion. We use the term interpretable, as opposed to inductive, to highlight that our approach, like sparse linear regression, gives actionable insights on what are the important features of matrix A. We note that Interpretable Matrix Completion is considerably harder than inductive or the classical matrix completion problem as it is a discrete problem and selecting k out of p factors is exponential in complexity.

In this paper, we show that the Interpretable Matrix Completion problem can be written as a mixed integer convex optimization problem. Inspired by Bertsimas and van Parys (2020) for sparse linear regression, we reformulate the interpretable matrix completion problem as a binary convex optimization problem. Then we introduce a new algorithm OptComplete, based on stochastic cutting planes, to enable scalability for matrices of sizes on the order of (). In addition, we provide empirical evidence on both synthetic and real-world data that OptComplete is able to match or exceed current state-of-the-art methods for inductive and general matrix completion on both speed and accuracy, despite OptComplete solving a more difficult problem.

Specifically, our contributions in this paper are as follows:

1. We introduce the interpretable matrix completion problem, and reformulate it as a binary

convex optimization problem that can be solved using cutting planes methods.

2. We propose a new novel approach to cutting planes by introducing stochastic cutting planes.

We prove that the new algorithm converges to an optimal solution of the interpretable matrix completion problem with exponentially vanishing failure probability.

3. We present computational results on both synthetic and real datasets that show that the

algorithm matches or outperforms current state-of-the-art methods in terms of both scalability and accuracy.

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 5

The structure of the paper is as follows. In Section 2, we introduce the binary convex reformulation of the low-rank interpretable matrix completion problem, and how it can be solved through a cutting plane algorithm, which we denote CutPlanes. In Section 3, we introduce OptComplete, a stochastic cutting planes method designed to scale the CutPlanes algorithm in Section 2, and show that it recovers the optimal solution of CutPlanes with exponentially vanishing failure probability. In Section 4, we report on computational experiments with synthetic data that compare OptComplete to Inductive Matrix Completion (IMC) introduced in Natarajan and Dhillon (2014) and SoftImpute-ALS (SIALS) by Hastie et al. (2015), two state-of-the-art matrix completion algorithms for inductive and general completion. We also compare OptComplete to CutPlanes to demonstrate the 20x to 60x speedup of the stochastic algorithm. In Section 5, we report computational experiments on the Netflix Prize dataset. In Section 6 we provide our conclusions.

Matrix completion has been applied successfully to many tasks, including recommender systems Koren et al. (2009), social network analysis Chiang et al. (2014) and clustering Chen et al. (2014b). After Cand`es and Tao (2010) proved a theoretical guarantee for the retrieval of the exact matrix under the nuclear norm convex relaxation, a lot of methods have focused on the nuclear norm problem (see Mazumder et al. (2010), Beck and Teboulle (2009), Jain et al. (2010), and Tanner and Wei (2013) for examples). Alternative methods include alternating projections by Recht and R´e (2013) and Grassmann manifold optimization by Keshavan et al. (2009). There has also been work where the uniform distributional assumptions required by the theoretical guarantees are violated, such as Negahban and Wainwright (2012) and Chen et al. (2014a).

Interest in inductive matrix completion intensified after Xu et al. (2013) showed that given predictive side information, one only needs O(logn) samples to retrieve the full matrix. Thus, most of this work (see Xu et al. (2013), Jain and Dhillon (2013), Farhat et al. (2013),

6 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

Natarajan and Dhillon (2014)) have focused on the case in which the side information is assumed to be perfectly predictive so that the theoretical bound of O(logn) sample complexity Xu et al. (2013) can be achieved. Chiang et al. (2015) explored the case in which the side information is corrupted with noise, while Shah et al. (2017) and Si et al. (2016) incorporated nonlinear combination of factors into the side information. Surprisingly, as pointed out by a recent article Nazarov et al. (2018), there is a considerable lack of effort to introduce sparsity/interpretability into inductive matrix completion, with Lu et al. (2016), Soni et al. (2016) and Nazarov et al. (2018) being among the only works that attempt to do so. Our work differs from the previous attempts in that previous attempts mainly focus on choosing latent features which are sparse linear combinations of the given features. In contrast interpretable matrix completion is aimed to select exactly k features from the known features.

2. Interpretable Matrix Completion

In this section, we present the mathematical formulation of Interpretable Matrix Completion and how it can be reformulated as a binary convex problem that is based on Bertsimas and van Parys (2020). We show how this naturally leads to a cutting plane algorithm, and discuss its computational complexity. We also discuss the two-sided information case, and how that reduces to the sparse regression problem.

2.1. Binary Convex Reformulation of Interpretable Matrix Completion

The (one-sided) interpretable matrix completion problem can be written as a mixed binary optimization problem:

where and:

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 7

We note that given that , the rank of matrix X is indeed k. We further note that the coefficients of S can be taken to be binary without loss of generality, since if they are not and ), then by applying the transformation:

for ), results in an equivalent problem with the coefficients of S being binary.

For this paper, we consider the squared norm, and for robustness purposes (see Bertsimas and van Parys (2020) and Bertsimas and Copenhaver (2018)), we add a Tikhonov regularization term to the original problem. Specifically, the (one-sided) interpretable matrix completion

problem with regularization we address is

In this section, we show how that problem (4) can be reformulated as a binary convex optimization problem, and can be solved to optimality using a cutting plane algorithm. The main theorem and proof is presented below:

, where is the ith row of A with unknown entries taken to be 0, and with the jth column of B.

Proof: With the diagonal matrices defined above, we can rewrite the sum in (4) over known entries of A, , as a sum over the rows of A:

8 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

where is the ith row of X. Using , then where is

the ith row of U. Moreover,

Then, Problem (4) becomes:

We then notice that within the sum each row of U can be optimized separately, leading to:

The inner optimization problem min+ 1can be solved in closed form given S, as it is a weighted linear regression problem with Tiknorov regularization, see Bertsimas and van Parys (2020). The closed form solution is:

So Problem (5) can be simplified to:

Finally, we notice that

and we obtain the required expression. Since are positive semi-definite, and the inverse of positive semi-definite matrices is a convex function, the entire function is convex in

With Theorem 1, our original problem can now be restated as:

This can be solved utilizing the cutting plane algorithm first introduced by Duran and Grossmann (1986), summarized as Algorithm 1.

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 9

The cutting plane algorithm, at iteration t, adds a linear approximation of c(s) at the current feasible solution to the set of constraints:

and we solve the mixed-integer linear programming problem:

to obtain . We see that is exactly the minimum value of the current approximation of ), defined below:

10 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

Since c(s) is convex, the piecewise linear approximation ) is an outer approximation (), so . As the algorithm progresses, the set of linear approximations form an increasingly better approximation of c(s), and increases with t. The algorithm terminates once does not further increase, as it implies the linear approximation shares the same minimum value as the true function c(s), which is the desired value.

Once the optimal solution is reached, we can obtain the optimal U using the closed form solution in (6) and recover X. In the next section, we discuss how this algorithm can be implemented in the context of c(s) in (7) and derive its computational complexity.

2.2. Implementation and Computational Complexity of CutPlanes

The computational complexity of the cutting plane comes from calculating c(s) and its derivative ). We first introduce the notations and .

Then, the function c(s) in (7) can be expressed as

To calculate the derivative ), it is easier to utilize the expression in Theorem 1 and then utilize

the chain rule. After some manipulations, we obtain

Therefore, we would focus on calculating ). First, by the Matrix Inversion Lemma (Woodbury

(1949)) we have

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 11

where is the feature matrix formed by the k columns of B such that have suppressed the dependency of V on s for notation ease. Note that in order to compute ) using Eq. (9) we need to invert an matrix, while from Eq. (12) we need to invert a matrix + , which only requires ) calculations. Furthermore, note that calculating

we do not need to multiply on the unknown entries. Similarly, we can compute

multiplications, and in multiplications.

Therefore, we can compute ) in floating point complexity of + ). Then to cal-

culate

respectively. Thus, the total complexity of generating a full cutting plane is:

2.3. Two-sided Information Case

In this section, we briefly discuss the matrix completion problem under the two-sided information case, and how it reduces to the problem of sparse linear regression. The two sided interpretable matrix completion problem with Tikhonov regularization can be stated as follows:

where is a known matrix of features of each row, is a known matrix of features of each column, and is a sparse matrix that has k nonzero entries, ensuring that Rank(. We note that in Eq. (14) we restrict the support of matrix L to be k, rather than forcing the entries of L to be binary. This is because unlike in the one-sided case, both U and B are known, so we cannot apply the scaling transformation in (3).

We denote by the ith column of U and the jth column of B. We introduce

the matrices as in Theorem 1. Using X = ULB, we can write

12 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

where is the (i,j)th entry of the matrix formed by multiplying qth column of U with th column of B. Then, Problem (14) becomes:

As every D matrix is known, this becomes a sparse regression problem where there are features to choose from (the D matrices), there are samples (the A matrix), the sparsity requirement is k, the regression coefficients are L, and we have Tikhonov regularization. Vectorizing D, L, and A reduces the problem back to the familiar form of sparse linear regression, that can be solved by the algorithm developed in Bertsimas and van Parys (2020) at scale.

3. OptComplete: The Stochastic Cutting Plane Speedup

In this section, we introduce OptComplete, a stochastic version of the cutting plane algorithm introduced in Section 2. We present theoretical results to show that the stochastic algorithm recovers the true optimal solution of the original algorithm with high probability without distributional assumptions. We also include a discussion on the dependence of such probability with various factors and its favorable theoretical computational complexity.

3.1. Introduction of OptComplete

In the previous section, we showed that through careful evaluation, we can calculate a full cutting plane in ) + ) calculations. However, in very high dimensions where are extremely large, the cost of generating the full cutting plane is still prohibitive. Thus, we consider generating approximations of the cutting plane that would enable the algorithm to scale for high values for n and m. Specifically, consider the cutting plane function in (10), reproduced below:

where:

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 13

We approximate the inner term ) by choosing 1 samples from {1,...,m} without replacement, with the set denoted F. Then we formulate the submatrix with such selected rows,

and similarly with

Then we choose 1 samples from {1,...,n} without replacement, with the set denoted G. We can then calculate an approximation of c(s) using the approximated (s):

where the set F is chosen independently for every row . Then the derivative of ˜) is:

Using such approximations, we can derive a stochastic cutting-plane algorithm, which we call OptComplete presented as Algorithm 2.

For this algorithm to work, we need the approximation ˜) and its derivative to be close to the nominal values. Furthermore, the approximated cutting planes should not cutoff the true solution. In the next section, we show that OptComplete enjoys such properties with high probability. In Section 3.3, we discuss how to select the size of f and g.

3.2. Main Theoretical Results

We would first show that the inner approximation is close to the true term with high probability:

Theorem 2 Let A be a partially known matrix, B a known feature matrix, and as defined in

Theorem 1. Let F be a random sample of size f from the set {1,...,m}, chosen without replacement.

14 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

We see that, without assumptions on the data, the inner approximation for both the value and

the derivative follows a bound with

inverting the statements give that, for all and all :

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 15

So the failure probability drops off exponentially with increasing bound , reflecting a Gaussian tail structure for ) and (s). The proof is contained in Appendix A.

Using this result, we are able to prove a tight deviation bound for the approximated cost function ) and ):

Theorem 3 Let A be a partially known matrix, B a known feature matrix, and as defined in Theorem 1. Let G be a random sample of size g from {1,...,n} chosen without replacement. Then for each , we let be a random sample of size f from the set {1,...,m}, all chosen without replacement. We have, with probability at least 1 :

Similar to the inner approximations, ) and ) has Gaussian tails. Furthermore, the scaling here only depends on g and not f: This shows that the error of the inner approximation is dominated by the outer sampling of the rows G. The proof is contained in Appendix B.

Then, using this result, we are able to prove our main result for OptComplete. We would first introduce a new definition:

Definition 1 The convexity parameter a of the cost function c(s) is defined as the largest positive number for which the the following statement is true:

We have the following proposition which shows that unless the cost function is degenerate (i.e. different sets of k features doesn’t change the solution), we always have a positive convexity parameter:

16 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

The proof is contained in Appendix C. Now, we state our main theorem for OptComplete.

Theorem 4 For the matrix completion problem (4), let be a known feature matrix, a matrix with entries partially known, and OptComplete as defined in Algorithm 2. Assume that Problem (4) is feasible. Then, OptComplete terminates in a finite number of steps C, and finds an optimal solution of (4) with probability at least 1 expwhere D is an absolute constant independent of C,f,g,k,p, and a is the convexity parameter of the functions

The proof is contained in Appendix D. This theorem shows that as long as the original problem is feasible, OptComplete is able to find the optimal solution of the original binary convex problem with exponentially vanishing failure probability that scales as exp. The theorem requires no assumptions on the data, and thus applies generally. We again note that the bound does not depend on f and only on g: we would discuss how this would inform our selection of the size of f and g in the next section.

3.3. Sampling Size and Computational Complexity

To select an appropriate f and g, we first note that Cand`es and Tao (2010) showed that to complete a square matrix of rank k, we need at least O(kN logN) elements. Assume an average known rate of in the original matrix A, the expected number of known elements under a sampling of f and g is . Using , we need that:

for some constant c. Theorem 4 showed that the bound on failure probability scales with exp, and thus we cannot have g too small. Using (13), the complexity of the cutting plane with f and g samples are:

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 17

Therefore, if we fix the expected known elements () constant, it is more advantageous to select a smaller g, as g scales with . Thus, we set:

Experimentally, we found

2). Therefore, by (18), the approximated cutting plane has a computational complexity of:

This scales in a square root fashion in n and m, rather than linearly in n and m for the full cutting plane. This allows OptComplete to enjoy a considerable speedup compared to CutPlanes, as demonstrated in Section 4.

4. Synthetic Data Experiments

We assume that the matrix A = UV + E, where , and E is an error matrix with individual elements sampled from N(0,0.01). We sample the elements of U and V from a uniform distribution of [0,1], and then randomly select a fraction to be missing. We formulate the feature matrix B by combining with a confounding matrix that contains unnecessary factors sampled similarly from the Uniform [0,1] distribution. We run OptComplete on a server with 16 CPU cores, using Gurobi 8.1.0. For each combination (), we ran 10 tests and report the median value for every statistic.

We report the following statistics with being the ground-truth factor vector, and

mated factor vector.

• n,m - the dimensions of A.

• p - the number of features in the feature matrix.

• k - the true number of features.

• - The fraction of missing entries in A.

• T - the total time taken for the algorithm.

18 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

• MAPE - the Mean Absolute Percentage Error (MAPE) for the retrieved matrix ˆA:

where is the set of missing data in A.

Since the concept of Interpretable Matrix Completion is new, there is a lack of directly comparable algorithms in the literature. Thus, in lieu, we compare OptComplete to state-of-the-art solvers for Inductive Matrix Completion and general matrix completion, which are:

• IMC by Natarajan and Dhillon (2014) - This algorithm is a well-accepted benchmark for

testing Inductive Matrix Completion algorithms.

• SoftImpute-ALS (SIALS) by Hastie et al. (2015) - This is widely recognized as a state-of-the-

art matrix completion method without feature information. It has among the best scaling behavior across all classes of matrix completion algorithms as it utilizes fast alternating least squares to achieve scalability.

We use the best existing implementations of IMC (Matlab 2018b) and SIALS (R 3.4.4, package softImpute) with parallelization on the same server.

We further compare our algorithm to CutPlanes, the original cutting plane algorithm developed in Section 2. It is known that for general mixed-integer convex problems, the cutting plane algorithm has the best overall performance (see e.g. Lubin et al. (2016) for details), and thus CutPlanes represent a good baseline of comparison for OptComplete.

We randomly selected 20% of those elements masked to serve as a validation set. The regularization parameter of OptComplete, the rank parameter of IMC and the penalization parameter of IMC and SIALS are selected using the validation set. The results are separated into sections below. The first five sections modify one single variable out of to investigate OptComplete’s scalability, where the leftmost column indicates the variable modified. The last section compares the four algorithms scalability for a variety of parameters that reflect more realistic scenarios.

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 19

Table 1 Comparison of OptComplete, IMC and SIALS on synthetic data. N/A means the algorithm did not

Overall, we see that OptComplete achieves near-exact retrieval on all datasets evaluated, and successfully recovers the factors in the ground truth. The solutions (and its error) also matches with that of CutPlanes, the standard cutting plane algorithm. The non-zero MAPE is due to the random noise added resulting in slightly perturbed coefficients.

20 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

For the realistic and large data sizes in the last panel, we see that OptComplete not only achieves near-exact retrieval, it does so while requiring considerably less time than IMC and SIALS at the same time. For m,n on the scale of and , OptComplete is over 20 times faster than IMC and over 40 times faster than SIALS. At the scale of and , IMC and SIALS did not finish running within 20 hours, while OptComplete completed in just under 12 minutes. We also see that OptComplete achieves very significant speedups compared to the standard cutting plane algorithm - up to 60x at the scale of and .

We analyze the scaling of OptComplete as a function of:

1. - The algorithm is able to retrieve the exact factors used even with 95% of missing data.

Furthermore, the running time decreased with increasing missing entries, consistent with the fact that is computational complexity scales with .

2. n - The algorithm has good scalability in n, reflecting its log(n)) type complexity. This

allows the algorithm to support matrices with n in the 10range. Its scaling behavior is superior to both IMC and SIALS.

3. m - The algorithm scales exceptionally well in m. We observe that empirically the algorithm

runtime seems to grow much slower than the theoretical log(m)) dependence. A closer examination reveals that as m increases, the number of cutting planes generated by Gurobi is decreasing. Qualitatively, this can be explained by a larger m giving the algorithm more signal to find which k features are the correct ones out of the p ones. We note that such behavior is also exhibited by CutPlanes, as it roughly scales as ) rather than the O(m) as expected. We see that IMC and SIALS scales as O(m).

4. p - The algorithm scales relatively well in p, which reflects the performance of the Gurobi

solver. We empirically observe that Gurobi is generating roughly O(1)) cutting planes. Thus, as each cutting plane is O(p), we expect ) dependence, as is observed here. We note that OptComplete achieves similar scaling behavior as IMC in p. Note here the SIALS algorithm does not utilize feature information and thus a change in p does not affect the algorithm’s run speed.

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 21

5. k - The algorithm does not scale very well in k. We empirically observe that Gurobi solver

is roughly generating O(k) cutting planes and each cutting plane has cubic dependence on k. It appears that SIALS and IMC almost have a linear scaling behavior. However, in most applications, such as recommendation systems or low-rank retrieval, k is usually kept very low (30), so this is not a particular concern.

5. Real-World Experiments

In this section, we report on the performance of OptComplete on the Netflix Prize dataset (Bennett et al. 2007). This dataset was released in a competition to predict ratings of customers on unseen movies, given over 10 million ratings scattered across 500,000 people and 16,000 movies. Thus, when presented in a matrix A where represents the rating of individual i on movie j, the goal is to complete the matrix A under a low-rank assumption.

The feature matrix B of OptComplete is constructed using data from the TMDB Database, and covers 59 features that measure geography, popularity, top actors/actresses, box office, runtime, genre and more. The full list of 59 features is contained in Appendix E.

For this experiment, we included movies where all 59 features are available, and people who had at least 5 ratings present. This gives a matrix of 471,268 people and 14,538 movies. The slight reduction of size from the original data is due to the lack of features for about 2,000 niche movies. To observe the scalability of OptComplete, we created five data sets:

1. Base - has dimensions 3,923 103.

2. Small - has dimensions 18,227 323.

3. Medium - has dimensions 96,601 788.

4. Large - has dimensions 471,268 1760.

5. Full - A has dimensions 471,268 14,538.

These sizes are constructed such that the total number of elements in A in the successive sizes are approximately different by approximately an order of magnitude.

22 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

For each individual matrix, we uniformly randomly withhold 20% of the ratings as a test set S, and use the remaining 80% of ratings to impute a complete matrix ˆA - we perform cross-validation on the appropriate hyperparameters. Then, we report MAPE.

For comparison, we again use IMC and SIALS. We set the maximum rank of SIALS to be k -the rank optimized for in OptComplete. The results are listed below:

Table 2 Comparison of methods on Netflix data for k = 5.

Table 3 Comparison of methods on Netflix data for k = 10.

We can see that OptComplete outperforms both IMC and SIALS in accuracy across the datasets under different k; furthermore in the two largest datasets OptComplete ran 10x to 20x faster than IMC and SIALS. Here we see that an increase from k = 5 to k = 10 actually decreased out-of-sample performance as additional factors are actually not very helpful in predictive customer tastes. The decline for OptComplete and IMC were especially higher due to the fact that the possible factors are fixed and thus an increase in the number of factors caused some non-predictive factors to be included.

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 23

For the k = 5 case, OptComplete identified the following as the top factors that influences an individual’s rating:

• IMDB Rating

• Genre: Drama

• Released within last 10 years

• Number of Top 100 Actors

• Produced in US

These factors provide an intuitive explanation of the individual ratings of each customer in terms of a small number of factors, while exceeding the high predictive accuracy of SIALS.

6. Conclusions

We have presented OptComplete, a scalable algorithm to retrieve a low-rank matrix in the presence of side information. Compared with state of the art algorithms for matrix completion, OptComplete exceeds current benchmarks on both scalability and accuracy and provides insight on the factors that affect the ratings.

24 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

Appendices A. Proof of Theorem 2

We first note that since is a finite set, we only need to prove the result for a particular , and it would apply for all s. Therefore, we would assume s is fixed below. For simplicity, we would

only demonstrate the proof for

exact fashion. Furthermore, since we are only focusing on one particular , we would drop all i subscripts below for ease of notation. The quantities of interest are therefore:

First, let

has orthogonal columns such that , and . Note such definition implies (1). Then, we would rewrite the terms as follows:

We note that

Therefore, if we treat

drawn without replacement from that set, and similarly for

Hoeffding’s inequality to bound the deviation of these terms, as reproduced below:

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 25

For a proof, see for example Boucheron et al. (2013). Then, applying Proposition 2 to

and inverting the inequality, we have

where A,B are constants independent of .

Now we would show that is close to:

Lemma 1

To prove this, we would first introduce a matrix analog of the well-known Chernoff bound, the proof of which can be found in Tropp (2012):

where is the minimum/maximum eigenvalue function. Sample uniformly at random without replacement. Compute:

26 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

Using this lemma, we would proceed with the proof of Lemma 1.

Proof of Lemma 1: First, we note that

where rank-one positive semi-definite matrices. Therefore, we can take as a random sample of size f from the set , which satisfies the conditions in Lemma

2 with D = O(k). Furthermore, with X, we observe that we have , so we have

Therefore, we apply Lemma 2 to and obtain

where we set with (1). Some rearrangement gives:

Now since , we have

Combining equation (A6) and (A5) gives

Then, we have

Taking log(2) gives the required result. (Note (1) as we setup the QR decomposition to have )).

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 27

With equations (A2), (A3), and (A4), we are now ready to bound ) and ). We first introduce another lemma from matrix perturbation theory (for proof, see e.g. Stewart (1990)).

Then, we have

Using (A2) and triangle inequality, we have

Using (A3) and Lemma 3, we have

for some O(1) constant . Note that , so using (A4), we obtain

for some new O(1) constant . Then, using (A10), we have

28 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

Since 1, we have

For some O(1) constant M with probability at least 1 . Therefore, taking and absorbing the extra constant into M gives the requires result.

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 29

B. Proof of Theorem 3

We first prove a Hoeffding-type bound as follows

Proposition 3 Let be independent (but not necessarily identically distributed) random

We note that such statement differs from Hoeffding’s inequality as we do not require . The proof is as follows.

Proof of Proposition 3: We first introduce a lemma known as Chernoff’s bounding method (proven in Chernoff (1952)):

Let . Then, we have

We aim to bound ], subject to exp

of ] is maximized when

function of exp

30 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

where erf is the standard error function. As the error function is upper bounded by 1, the last expression is less than or equal to:

We then substitute this result into (A12) and obtain

Note that this is minimized at , so we have

Therefore,

By applying the previous derivation to , we obtain

Combining the two results completes the proof.

Using Proposition 3, we bound the difference between ˜) and ). We have

The first term can be seen as the tail bound for a random sample of size g chosen without replacement from the finite set . Thus, we can apply Hoeffding’s theorem in Proposition 2,

and obtain that with probability at least 1

Substituting (A13) into the expression above shows that with probability at least 1 we have:

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 31

Note that, for any fixed set G, for with , we have that (s) and (s) are independent

(as we construct F separately for every i). Furthermore, Theorem 2 can be inverted to read

). Then, applying Proposition 3 to the second term, we have

with probability 1 . As 1, a loose bound is therefore

with probability 1 . Therefore, taking and log(2) we have

Through a similar derivation, we have

32 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

C. Proof of Proposition 1

We would prove the contrapositive. Assume that a = 0. We would calculate the second derivative of the cost function expression in 7. First, define + . Then our cost function can be rewritten as:

Then the second derivative can be easily calculated as:

Thus, for any vector , we have:

Now since a = 0, the Hessian cannot be positive definite for all s (as if it was, then a > 0). Then there exist such that for , we have:

and t Note that for any 0, is positive definite, so is positive definite for any

for all i. Now note that we have:

We have that:

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 33

And similarly, we have

Thus, we have:

Therefore, if we have a = 0, then we must have equality on the cost function for two feasible solutions. Thus, if the cost function is not degenerate for feasible solutions, then a > 0.

34 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

D. Proof of Theorem 4

In this section, we provide the proof of Theorem 4. We first note that OptComplete is a specific implementation of the outer approximation algorithm, which Fletcher and Leyffer (1994) have shown to always terminate in finite number of steps C. Thus, given that we assumed the problem is feasible, for OptComplete to not return an optimal solution, it would have to cut it off during the course of its execution. Let be an optimal solution for Problem (4). Let be an optimal solution at the t-th iteration of OptComplete, ]. The cutting plane constraint for OptComplete at the

point of an optimal solution is

If , then will be cut off, and OptComplete will not find . Applying the definition of

the convexity parameter (16) and letting (noting that 1) we obtain

Therefore, if

or equivalently if

then OptComplete will not find . Therefore, for OptComplete to succeed, should satisfy 2 for all

Let . Then,

Since at each step of OptComplete we randomly sample r new rows and s new columns, the events

are independent for different ], and hence

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 35

Therefore, we focus on calculating 2). Since

using Theorem 3, we can thus provide the following bound for deviation of , for some constant

where . Then we can invert this to calculate 2)

completing the proof.

36 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

E. List of Features Used in the Netﬂix Problem

• 24 Indicator Variables for Genres: Action, Adventure, Animation, Biography, Comedy, Crime,

Documentary, Drama, Family, Fantasy, Film Noir, History, Horror, Music, Musical, Mystery,

Romance, Sci-Fi, Short, Sport, Superhero, Thriller, War, Western

• 5 Indicator Variables for Release Date: Within last 10 years, Between 10-20 years, Between

20-30 years, Between 30-40 years, Between 40-50 Years

• 6 Indicator Variables for Top Actors/Actresses defined by their Influence Score at time of

release: Top 100 Actors, Top 100 Actresses, Top 250 Actors, Top 250 Actresses, Top 1000 Actors,

Top 1000 Actresses

• IMDB Rating

• Number of Reviews

• Total Production Budget

• Total Runtime

• Total Box Office Revenue

• Indicator Variable for whether it is US produced

• 11 Indicator Variables for Month of Year Released (January removed to prevent multicollinear-

ity)

• Number of Original Music Score

• Number of Male Actors

• Number of Female Factors

• 3 Indicator Variables for Film Language: English, French, Japanese

• Constant

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 37

References

Beck A, Teboulle M (2009) A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences 2(1):183–202.

Bennett J, Lanning S, et al. (2007) The netflix prize. KDD Cup and Workshop 2007 (Citeseer).

Bertsimas D, Copenhaver MS (2018) Characterization of the equivalence of robustification and regularization in linear, median, and matrix regression. European Journal of Operations Research 270:931–942.

Bertsimas D, van Parys B (2020) Sparse high-dimensional regression: Exact scalable algorithms and phase transitions. Annals of Statistics 48(1):300–323.

Boucheron S, Lugosi G, Massart P (2013) Concentration inequalities: A nonasymptotic theory of independence (Oxford university press).

Boyd S, El Ghaoui L, Feron E, Balakrishnan V (1994) Linear matrix inequalities in system and control theory, volume 15 (SIAM).

Candes EJ, Plan Y (2010) Matrix completion with noise. Proceedings of the IEEE 98(6):925–936.

Cand`es EJ, Tao T (2010) The power of convex relaxation: Near-optimal matrix completion. IEEE Transactions on Information Theory 56(5):2053–2080.

Chen Y, Bhojanapalli S, Sanghavi S, Ward R (2014a) Coherent matrix completion. International Conference on Machine Learning, 674–682.

Chen Y, Jalali A, Sanghavi S, Xu H (2014b) Clustering partially observed graphs via convex optimization. The Journal of Machine Learning Research 15(1):2213–2238.

Chernoff H (1952) A measure of asymptotic efficiency for tests of a hypothesis based on the sum of observa- tions. The Annals of Mathematical Statistics 23(4):493–507.

Chiang KY, Hsieh CJ, Dhillon IS (2015) Matrix completion with noisy side information. Advances in Neural Information Processing Systems, 3447–3455.

Chiang KY, Hsieh CJ, Natarajan N, Dhillon IS, Tewari A (2014) Prediction and clustering in signed networks: a local to global perspective. The Journal of Machine Learning Research 15(1):1177–1213.

Dhillon P, Lu Y, Foster DP, Ungar L (2013) New subsampling algorithms for fast least squares regression. Advances in Neural Information Processing Systems, 360–368.

38 Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

Duran MA, Grossmann IE (1986) An outer-approximation algorithm for a class of mixed-integer nonlinear programs. Mathematical Programming 36(3):307–339.

Farhat MR, Shapiro BJ, Kieser KJ, Sultana R, Jacobson KR, Victor TC, Warren RM, Streicher EM, Calver A, Sloutsky A, et al. (2013) Genomic analysis identifies targets of convergent positive selection in drug-resistant mycobacterium tuberculosis. Nature Genetics 45(10):1183.

Fletcher R, Leyffer S (1994) Solving mixed integer nonlinear programs by outer approximation. Mathematical programming 66(1-3):327–349.

Hastie T, Mazumder R, Lee JD, Zadeh R (2015) Matrix completion and low-rank svd via fast alternating least squares. The Journal of Machine Learning Research 16(1):3367–3402.

Jain P, Dhillon IS (2013) Provable inductive matrix completion. arXiv preprint arXiv:1306.0626 .

Jain P, Meka R, Dhillon IS (2010) Guaranteed rank minimization via singular value projection. Advances in Neural Information Processing Systems, 937–945.

Ji H, Liu C, Shen Z, Xu Y (2010) Robust video denoising using low rank matrix completion. Computer Society Conference on Computer Vision and Pattern Recognition (IEEE).

Keshavan RH, Oh S, Montanari A (2009) Matrix completion from a few entries. Information Theory, 2009. ISIT 2009. IEEE International Symposium on, 324–328 (IEEE).

Koren Y, Bell R, Volinsky C (2009) Matrix factorization techniques for recommender systems. Computer 30–37.

Lu J, Liang G, Sun J, Bi J (2016) A sparse interactive model for matrix completion with side information. Advances in Neural Information Processing Ssystems, 4071–4079.

Lubin M, Yamangil E, Bent R, Vielma JP (2016) Extended formulations in mixed-integer convex program- ming. International Conference on Integer Programming and Combinatorial Optimization, 102–113 (Springer).

Mazumder R, Hastie T, Tibshirani R (2010) Spectral regularization algorithms for learning large incomplete matrices. Journal of Machine Learning Research 11(Aug):2287–2322.

Natarajan N, Dhillon IS (2014) Inductive matrix completion for predicting gene–disease associations. Bioinformatics 30(12):160–168.

Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!) 39

Nazarov I, Shirokikh B, Burkina M, Fedonin G, Panov M (2018) Sparse group inductive matrix completion. arXiv preprint arXiv:1804.10653 .

Negahban S, Wainwright MJ (2012) Restricted strong convexity and weighted matrix completion: Optimal bounds with noise. Journal of Machine Learning Research 13(May):1665–1697.

Recht B, R´e C (2013) Parallel stochastic gradient algorithms for large-scale matrix completion. Mathematical Programming Computation 5(2):201–226.

Shah V, Rao N, Ding W (2017) Matrix factorization with side and higher order information. Stat 1050:4.

Si S, Chiang KY, Hsieh CJ, Rao N, Dhillon IS (2016) Goal-directed inductive matrix completion. Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 1165– 1174 (ACM).

Soni A, Chevalier T, Jain S (2016) Noisy inductive matrix completion under sparse factor models. arXiv preprint arXiv:1609.03958 .

Stewart GW (1990) Matrix perturbation theory .

Tanner J, Wei K (2013) Normalized iterative hard thresholding for matrix completion. SIAM Journal on Scientific Computing 35(5):S104–S125.

Tropp JA (2012) User-friendly tail bounds for sums of random matrices. Foundations of computational mathematics 12(4):389–434.

Woodbury MA (1949) The stability of out-input matrices. Chicago, Ill .

Xu M, Jin R, Zhou ZH (2013) Speedup matrix completion with side information: Application to multi-label learning. Advances in Neural Information Processing Systems, 2301–2309.

Designed for Accessibility and to further Open Science