Interpretable Matrix Completion: A Discrete Optimization Approach

TL;DR

This paper introduces Interpretable Matrix Completion, a novel approach that uses side information and a binary convex optimization reformulation to provide meaningful insights into low-rank matrices, scalable to very large sizes.

Contribution

It presents OptComplete, a scalable algorithm based on stochastic cutting planes for interpretable matrix completion, with demonstrated efficiency and accuracy on large datasets.

Findings

OptComplete scales efficiently to matrices of size 10^6 by 10^6.

The method achieves competitive accuracy compared to state-of-the-art techniques.

It provides interpretability by revealing factors influencing the matrix.

Abstract

We consider the problem of matrix completion on an $n \times m$ matrix. We introduce the problem of Interpretable Matrix Completion that aims to provide meaningful insights for the low-rank matrix using side information. We show that the problem can be reformulated as a binary convex optimization problem. We design OptComplete, based on a novel concept of stochastic cutting planes to enable efficient scaling of the algorithm up to matrices of sizes $n=10^6$ and $m=10^6$. We report experiments on both synthetic and real-world datasets that show that OptComplete has favorable scaling behavior and accuracy when compared with state-of-the-art methods for other types of matrix completion, while providing insight on the factors that affect the matrix.

Tables

Table 1. Table 1: Comparison of OptComplete, IMC and SIALS on synthetic data. N / A 𝑁 𝐴 N/A means the algorithm did not complete running in 20 hours, corresponding to 72000 seconds.

$𝒏$

$𝒎$

$𝒑$

$𝒌$

$𝝁\%$

OptComplete

CutPlanes

IMC

SIALS

$𝑻$

MAPE

$𝑻$

MAPE

$𝑻$

MAPE

$𝑻$

MAPE

$\mu$

100

100

15

5

$20\%$

1.7s

$0.1\%$

6.0s

$0.1\%$

0.03s

$0.01\%$

0.02s

$0.3\%$

100

100

15

5

$50\%$

0.9s

$0.02\%$

4.5s

$0.02\%$

0.07s

$0.5\%$

0.03s

$0.9\%$

100

100

15

5

$80\%$

0.6s

$0.03\%$

2.5s

$0.03\%$

0.09s

$1.3\%$

0.06s

$5.6\%$

100

100

15

5

$95\%$

0.2s

$0.04\%$

1.2s

$0.04\%$

0.12s

$12.1\%$

0.12s

$7.4\%$

$n$

100

100

15

5

$50\%$

0.9s

$0.02\%$

4.5s

$0.02\%$

0.07s

$0.5\%$

0.03s

$0.9\%$

${10}^3$

100

15

5

$50\%$

3.1s

$0.01\%$

72.5s

$0.01\%$

0.6s

$0.4\%$

0.1s

$0.2\%$

${10}^4$

100

15

5

$50\%$

9.5s

$0.004\%$

957s

$0.004\%$

4.5s

$0.3\%$

6.5s

$0.5\%$

${10}^5$

100

15

5

$50\%$

18.0s

$0.003\%$

10856s

$0.003\%$

32.7s

$0.1\%$

38s

$3.0\%$

$m$

100

100

15

5

$50\%$

0.9s

$0.02\%$

4.5s

$0.02\%$

0.07s

$0.5\%$

0.03s

$0.9\%$

100

${10}^3$

15

5

$50\%$

0.7s

$0.01\%$

18.6s

$0.01\%$

0.8s

$0.3\%$

0.1s

$0.5\%$

100

${10}^4$

15

5

$50\%$

1.2s

$0.004\%$

68.5s

$0.004\%$

6.2s

$0.2\%$

0.8s

$0.3\%$

100

${10}^5$

15

5

$50\%$

3.0s

$0.002\%$

259s

$0.002\%$

56.2s

$0.1\%$

12.7s

$0.8\%$

$p$

100

100

15

5

$50\%$

0.9s

$0.02\%$

4.5s

$0.02\%$

0.07s

$0.5\%$

0.03s

$0.9\%$

100

100

50

5

$50\%$

2.0s

$0.02\%$

18.0s

$0.02\%$

0.3s

$0.6\%$

0.03s

$0.9\%$

100

100

200

5

$50\%$

12.1s

$0.02\%$

95.9s

$0.02\%$

1.9s

$0.8\%$

0.03s

$0.9\%$

100

100

${10}^3$

5

$50\%$

90.3s

$0.02\%$

680s

$0.02\%$

10.4s

$1.0\%$

0.03s

$0.9\%$

$k$

100

100

50

5

$50\%$

2.0s

$0.02\%$

18.0s

$0.02\%$

0.3s

$0.5\%$

0.03s

$0.9\%$

100

100

50

10

$50\%$

20.7s

$0.06\%$

130s

0.06%

0.20s

$1.2\%$

0.1s

$0.8\%$

100

100

50

20

$50\%$

240s

$0.07\%$

1584s

$0.07\%$

0.35s

$2.1\%$

0.21s

$1.0\%$

100

100

50

30

$50\%$

980s

$0.09\%$

8461s

$0.09\%$

0.5s

$3.3\%$

0.43s

$2.8\%$

100

100

15

5

$95\%$

0.2s

$0.04\%$

1.2s

$0.04\%$

0.12s

$12.1\%$

0.12s

$7.4\%$

${10}^3$

${10}^3$

50

5

$95\%$

1.4s

$0.006\%$

3.5s

$0.006\%$

$4.6s$

$4.7\%$

$2.8s$

$12.5\%$

${10}^4$

${10}^3$

100

5

$95\%$

5.7s

$0.002\%$

35.2s

$0.002\%$

$18s$

$2.5\%$

$20.7s$

$12.6\%$

${10}^5$

${10}^3$

200

10

$95\%$

52s

$0.001\%$

1520s

$0.001\%$

295s

$1.7\%$

$420s$

$4.6\%$

${10}^5$

${10}^4$

200

10

$95\%$

98s

$0.001\%$

5769s

$0.001\%$

$1750s$

$0.5\%$

$4042s$

$4.1\%$

${10}^6$

${10}^4$

200

10

$95\%$

480s

$0.001\%$

$N/A$

$N/A$

$13750s$

$0.3\%$

$25094s$

$2.5\%$

${10}^6$

${10}^5$

200

10

$95\%$

680s

$0.001\%$

$N/A$

$N/A$

$N/A$

$N/A$

$N/A$

$N/A$

${10}^6$

${10}^6$

200

10

$95\%$

1415s

$0.001\%$

$N/A$

$N/A$

$N/A$

$N/A$

$N/A$

$N/A$

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

$𝒏$	$𝒎$	$𝒑$	$𝒌$	$𝝁\%$	OptComplete		IMC		SIALS
					$𝑻$	MAPE	$𝑻$	MAPE	$𝑻$	MAPE
3,923	103	59	5	$92.6\%$	6.0s	$29.4\%$	0.6s	$34.2\%$	0.3s	$31.2\%$
18,227	323	59	5	$94.8\%$	12.2s	$21.8\%$	5.2s	$29.1\%$	4.1s	$24.1\%$
96,601	788	59	5	$94.2\%$	25.5s	$20.9\%$	38.1s	$28.7\%$	30.4s	$21.3\%$
471,268	1,760	59	5	$93.6\%$	102s	$18.8\%$	460s	$24.6\%$	430s	$19.8\%$
471,268	14,538	59	5	$94.1\%$	170s	$15.7\%$	3921s	$21.5\%$	5300s	$16.7\%$

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

$𝒏$	$𝒎$	$𝒑$	$𝒌$	$𝝁\%$	OptComplete		IMC		SIALS
					$𝑻$	MAPE	$𝑻$	MAPE	$𝑻$	MAPE
3,923	103	59	10	$92.6\%$	11.0s	$30.4\%$	1.4s	$36.7\%$	0.8s	$35.8\%$
18,227	323	59	10	$94.8\%$	20.3s	$24.0\%$	12.5s	$32.5\%$	7.0s	$28.9\%$
96,601	788	59	10	$94.2\%$	45.9s	$22.3\%$	84.2s	$29.6\%$	50.7s	$22.8\%$
471,268	1,760	59	10	$93.6\%$	260s	$20.7\%$	1022s	$24.8\%$	870s	$20.7\%$
471,268	14,538	59	10	$94.1\%$	380s	$19.6\%$	8704s	$23.1\%$	10240s	$20.0\%$

Equations

$$\underset{\mathbf{X}}{\min }\frac{1}{nm}\sum _{(i,j)\in \Omega }\Vert X_{ij}-A_{ij}\Vert \text{subject to}\mathrm{Rank}(\mathbf{X})=k,$$

$$\underset{\mathbf{U}}{\min }\underset{\mathbf{V}}{\min }\frac{1}{nm}\sum _{(i,j)\in \Omega }\Vert X_{ij}-A_{ij}\Vert \text{subject to}\mathbf{X}=\mathbf{U}{\mathbf{V}}^T.$$

$$\mathbf{V}=\mathbf{B}\mathbf{S}$$

$$\underset{\mathbf{U}}{\min }\underset{\mathbf{S}}{\min }\frac{1}{nm}\sum _{(i,j)\in \Omega }\Vert X_{ij}-A_{ij}\Vert \text{subject to}\mathbf{X}=\mathbf{U}{\mathbf{S}}^T{\mathbf{B}}^T.$$

$$5.6\times \text{running time}-0.00067\times \text{budget}+12\times \text{\# of Top 100 actors},$$

$$0.25\times \text{box office}-5\times \text{\# of Top 100 actors}.$$

$$\underset{\mathbf{U}}{\min }\underset{\mathbf{s}\in S_k^p}{\min }\frac{1}{nm}\sum _{(i,j)\in \Omega }\Vert X_{ij}-A_{ij}\Vert \text{subject to}\mathbf{X}=\mathbf{U}\mathbf{S}{\mathbf{B}}^T,$$

$$S_k^p=\left\{\mathbf{s}=(s_1,\dots ,s_p{)}^T\in \{0,1{\}}^p:\sum _{i=1}^p{s}_i=k\right\}.$$

$$\mathbf{U}\to \mathbf{U}\mathbf{D}\mathbf{S}\to \mathbf{S}{\mathbf{D}}^{-1}$$

$$\underset{\mathbf{U}}{\min }\underset{\mathbf{s}\in S_k^p}{\min }\frac{1}{nm}\left(\sum _{(i,j)\in \Omega }(X_{ij}-A_{ij}{)}^2+\frac{1}{\gamma }\Vert \mathbf{U}{\Vert }_2^2\right)\text{subject to}\mathbf{X}=\mathbf{U}\mathbf{S}{\mathbf{B}}^T.$$

$$\underset{\mathbf{s}\in S_k^p}{\min }c(\mathbf{s})=\frac{1}{nm}\sum _{i=1}^n{\overline{\mathbf{a}}}_i{\left({\mathbf{I}}_m+\gamma {\mathbf{W}}_i\left(\sum _{j=1}^p{s}_j{\mathbf{K}}_j\right){\mathbf{W}}_i\right)}^{-1}{\overline{\mathbf{a}}}_i^T,$$

$$({\mathbf{W}}_i{)}_{jj}=\{\begin{array}{ll}1,& (i,j)\in \Omega ,\text{vspace}3pt\\ 0,& (i,j)\notin \Omega ,\end{array}$$

$$\sum _{i=1}^n\Vert ({\mathbf{x}}_i-{\mathbf{a}}_i){\mathbf{W}}_i{\Vert }_2^2,$$

$$\Vert \mathbf{U}{\Vert }_2^2=\sum _{i=1}^n\Vert {\mathbf{u}}_i{\Vert }_2^2.$$

$$\underset{\mathbf{s}\in S_k^p}{\min }\underset{\mathbf{U}}{\min }\frac{1}{nm}\left(\sum _{i=1}^n\left(\Vert ({\mathbf{u}}_i\mathbf{S}{\mathbf{B}}^T-{\mathbf{a}}_i){\mathbf{W}}_i{\Vert }_2^2+\frac{1}{\gamma }\Vert {\mathbf{u}}_i{\Vert }_2^2\right)\right).$$

$$\underset{\mathbf{s}\in S_k^p}{\min }\frac{1}{nm}\left(\sum _{i=1}^n\underset{{\mathbf{u}}_i}{\min }\left(\Vert ({\mathbf{u}}_i\mathbf{S}{\mathbf{B}}^T-{\mathbf{a}}_i){\mathbf{W}}_i{\Vert }_2^2+\frac{1}{\gamma }\Vert {\mathbf{u}}_i{\Vert }_2^2\right)\right).$$

$$\underset{{\mathbf{u}}_i}{\min }\Vert ({\mathbf{u}}_i\mathbf{S}{\mathbf{B}}^T-{\mathbf{a}}_i){\mathbf{W}}_i{\Vert }_2^2+\frac{1}{\gamma }\Vert {\mathbf{u}}_i{\Vert }_2^2={\overline{\mathbf{a}}}_i({\mathbf{I}}_m+\gamma {\mathbf{W}}_i\mathbf{B}\mathbf{S}{\mathbf{B}}^T{\mathbf{W}}_i{)}^{-1}{\overline{\mathbf{a}}}_i^T.$$

$$\underset{\mathbf{s}\in S_k^p}{\min }\frac{1}{nm}\left(\sum _{i=1}^n{\overline{\mathbf{a}}}_i({\mathbf{I}}_m+\gamma {\mathbf{W}}_i\mathbf{B}\mathbf{S}{\mathbf{B}}^T{\mathbf{W}}_i{)}^{-1}{\overline{\mathbf{a}}}_i^T\right).$$

$$\mathbf{B}\mathbf{S}{\mathbf{B}}^T=\sum _{j=1}^p{s}_j{\mathbf{b}}^j({\mathbf{b}}^j{)}^T=\sum _{j=1}^p{s}_j{\mathbf{K}}_j,$$

$$\underset{\mathbf{s}\in S_k^p}{\min }c(\mathbf{s})=\frac{1}{nm}\sum _{i=1}^n{\overline{\mathbf{a}}}_i{\left({\mathbf{I}}_m+\gamma {\mathbf{W}}_i\left(\sum _{j=1}^p{s}_j{\mathbf{K}}_j\right){\mathbf{W}}_i\right)}^{-1}{\overline{\mathbf{a}}}_i^T.$$

$$\eta \ge c({\mathbf{s}}_t)+\nabla c({\mathbf{s}}_t{)}^T(\mathbf{s}-{\mathbf{s}}_t),$$

$$\underset{\mathbf{s}\in S_k^p,\eta \ge 0}{\min }\eta$$

$$\eta \ge c({\mathbf{s}}_t)+\nabla c({\mathbf{s}}_t{)}^T(\mathbf{s}-{\mathbf{s}}_t),i\in [t]$$

$${\eta }_{t+1}=\underset{\mathbf{s}}{\min }\underset{i\in [t]}{\max }c({\mathbf{s}}_t)+\nabla x({\mathbf{s}}_t{)}^T(\mathbf{s}-{\mathbf{s}}_t)=\underset{\mathbf{s}}{\min }{c}_t(\mathbf{s}).$$

$${\alpha }_i(\mathbf{s})=\frac{1}{m}{\overline{\mathbf{a}}}_i{\mathbf{\gamma }}_i(\mathbf{s})=\frac{1}{m}{\overline{\mathbf{a}}}_i\left[{\left({\mathbf{I}}_m+\gamma {\mathbf{W}}_i\left(\sum _{j=1}^p{s}_j{\mathbf{K}}_j\right){\mathbf{W}}_i\right)}^{-1}{\overline{\mathbf{a}}}_i^T\right],i=1,\dots ,n.$$

$$c(\mathbf{s})=\frac{1}{n}\sum _{i=1}^n{\alpha }_i(\mathbf{s})=\frac{1}{n}\sum _{i=1}^n\frac{{\overline{\mathbf{a}}}_i{\mathbf{\gamma }}_i(\mathbf{s})}{m}.$$

$$\nabla c(\mathbf{s})=\frac{1}{n}\sum _{i=1}^n-\frac{\gamma {\left({\mathbf{B}}^T{\mathbf{W}}_i{\mathbf{\gamma }}_i(\mathbf{s})\right)}^2}{m}.$$

$${\mathbf{\gamma }}_i(\mathbf{s})$$

$$=\left({\mathbf{I}}_m-\mathbf{V}{\left(\frac{{\mathbf{I}}_k}{\gamma }+{\mathbf{V}}^T{\mathbf{W}}_i\mathbf{V}\right)}^{-1}{\mathbf{V}}^T\right){\overline{\mathbf{a}}}_i^T$$

$$=\left({\overline{\mathbf{a}}}_i^T-\mathbf{V}{\left(\frac{{\mathbf{I}}_k}{\gamma }+{\mathbf{V}}^T{\mathbf{W}}_i\mathbf{V}\right)}^{-1}{\mathbf{V}}^T{\overline{\mathbf{a}}}_i^T\right),$$

$$\sum _{i=1}^{n}O(|{\Omega }_i|p+|{\Omega }_i|k^2+k^3)=O(|\Omega |(p+k^2)+n{k}^3).$$

$$\underset{\mathbf{L}}{\min }\frac{1}{nm}\left(\sum _{(i,j)\in \Omega }(X_{ij}-A_{ij}{)}^2+\frac{1}{\gamma }\Vert \mathbf{L}{\Vert }_2^2\right)\text{subject to}\mathbf{X}=\mathbf{U}\mathbf{L}{\mathbf{B}}^T\Vert \mathbf{L}{\Vert }_0=k,$$

Full text

\DoubleSpacedXI\TheoremsNumberedThrough\ECRepeatTheorems\EquationsNumberedThrough

\RUNAUTHOR

Bertsimas and Li

\RUNTITLE

Interpretable Matrix Completion

\TITLE

Interpretable Matrix Completion: A Discrete Optimization Approach

\ARTICLEAUTHORS\AUTHOR

Dimitris Bertsimas \AFFSloan School of Management and Operations Research Center, Massachusetts Institute of Technology, Cambridge, MA 02139, \EMAILdbertsim@mit.edu \AUTHORMichael Lingzhi Li \AFFOperations Research Center, Massachusetts Institute of Technology, Cambridge, MA 02139, \EMAILmlli@mit.edu

\ABSTRACT

We consider the problem of matrix completion on an $n\times m$ matrix. We introduce the problem of Interpretable Matrix Completion that aims to provide meaningful insights for the low-rank matrix using side information. We show that the problem can be reformulated as a binary convex optimization problem. We design OptComplete, based on a novel concept of stochastic cutting planes to enable efficient scaling of the algorithm up to matrices of sizes $n=10^{6}$ and $m=10^{6}$. We report experiments on both synthetic and real-world datasets that show that OptComplete has favorable scaling behavior and accuracy when compared with state-of-the-art methods for other types of matrix completion, while providing insight on the factors that affect the matrix.

\KEYWORDS

Matrix Completion, Mixed-Integer Optimization, Stochastic Approximation

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.

The classical low-rank matrix completion problem considers the following problem: Given a matrix $\bm{A}\in\mathbb{R}^{n\times m}$ with entries only partially known (denote $\Omega\subset\{1,\ldots,n\}\times\{1,\ldots,m\}$ as the set of known entries), we aim to recover a matrix $\bm{X}\in\mathbb{R}^{n\times m}$ of rank $k$ that minimizes a certain distance metric between $\bm{X}$ and $\bm{A}$ on the known entries of $\bm{A}$:

[TABLE]

where we normalized the objective so that it is $O(1)$. The rank $k$ constraint on $\bm{X}$ can be equivalently formulated as the existence of two matrices $\bm{U}\in\mathbb{R}^{n\times k}$, $\bm{V}\in\mathbb{R}^{m\times k}$ such that $\bm{X}=\bm{U}\bm{V}^{T}$. Therefore, the problem can be restated as:

[TABLE]

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 $A_{ij}$ being the response data of individual $i$ on item $j$. Therefore, the matrices $\bm{U}$ and $\bm{V}$ are commonly interpreted as the “user matrix” and “product matrix” respectively.

Let us denote each row of $\bm{U}$ as $\bm{u}_{i}$ and each column as $\bm{u}^{i}$ (similarly for $\bm{V}$). Then, $\bm{u}^{i}$ ($\bm{v}^{i}$) 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 $j$, $X_{ij}=\bm{u}_{i}\cdot\bm{v}_{j}$, 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 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 $\bm{B}\in\mathbb{R}^{m\times p}$. As an example, if the items are movies, then $\bm{b}^{j}$ represents feature $j$ for a movie (actors, budget, running time, etc), and each product feature $\bm{v}^{j}$ needs to be linear combination of such features. Mathematically, this translates to the constraint:

[TABLE]

where $\bm{S}\in\mathbb{R}^{p\times k}$. Therefore, the inductive version of the problem in (1) can be written as:

[TABLE]

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:

[TABLE]

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 significantly 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 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 $\bm{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 $(n,m)=(10^{6},10^{6})$. 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:

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. 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. 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.

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 $20$x to $60$x speedup of the stochastic algorithm. In Section 5, we report computational experiments on the Netflix Prize dataset. In Section 6 we provide our conclusions.

Literature

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ès 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é (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(\log n)$ samples to retrieve the full matrix. Thus, most of this work (see Xu et al. (2013), Jain and Dhillon (2013), Farhat et al. (2013), 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(\log n)$ 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:

[TABLE]

where $\bm{S}=\operatorname{Diag}\{s_{1},\ldots,s_{p}\}\in\mathbb{R}^{p\times p}$ and:

[TABLE]

We note that given that $\sum_{i=1}^{p}s_{i}=k$, the rank of matrix $\bm{X}$ is indeed $k$. We further note that the coefficients of $\bm{S}$ can be taken to be binary without loss of generality, since if they are not and $\bm{S}=\operatorname{Diag}(1/d_{1},\ldots,1/d_{p})$, then by applying the transformation:

[TABLE]

for $\bm{D}=\operatorname{Diag}(d_{1},\ldots,d_{p})$, results in an equivalent problem with the coefficients of $\bm{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

[TABLE]

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:

Theorem 1

Problem (4) can be reformulated as a binary convex optimization problem:

[TABLE]

where $\bm{W}_{1},\ldots,\bm{W}_{n}\in\mathbb{R}^{m\times m}$ are diagonal matrices:

[TABLE]

$\overline{\bm{a}}_{i}=\bm{W}_{i}\bm{a}_{i}$, $i=1\ldots,n$, where $\bm{a}_{i}\in\mathbb{R}^{1\times m}$ is the $i$th row of $\bm{A}$ with unknown entries taken to be 0, and $\bm{K}_{j}=\bm{b}^{j}(\bm{b}^{j})^{T}\in\mathbb{R}^{m\times m}$, $j=1,\ldots,p$ with $\bm{b}^{j}\in\mathbb{R}^{m\times 1}$ the $j$th column of $\bm{B}$.

Proof.

Proof: With the diagonal matrices $\bm{W}_{i}$ defined above, we can rewrite the sum in (4) over known entries of $\bm{A}$, $\sum_{(i,j)\in\Omega}(X_{ij}-A_{ij})^{2}$, as a sum over the rows of $\bm{A}$:

[TABLE]

where $\bm{x}_{i}\in\mathbb{R}^{1\times m}$ is the $i$th row of $\bm{X}$. Using $\bm{X}=\bm{U}\bm{S}\bm{B}^{T}$, then $\bm{x}_{i}=\bm{u}_{i}\bm{S}\bm{B}^{T}$ where $\bm{u}_{i}\in\mathbb{R}^{1\times m}$ is the $i$th row of $U$. Moreover,

[TABLE]

Then, Problem (4) becomes:

[TABLE]

We then notice that within the sum $\sum_{i=1}^{n}$ each row of $\bm{U}$ can be optimized separately, leading to:

[TABLE]

The inner optimization problem $\displaystyle\min_{\bm{u}_{i}}\|(\bm{u}_{i}\bm{S}\bm{B}^{T}-\bm{a}_{i})\bm{W}_{i}\|_{2}^{2}+\frac{1}{\gamma}\|\bm{u}_{i}\|_{2}^{2}$ can be solved in closed form given $\bm{S}$, as it is a weighted linear regression problem with Tiknorov regularization, see Bertsimas and van Parys (2020). The closed form solution is:

[TABLE]

So Problem (5) can be simplified to:

[TABLE]

Finally, we notice that

[TABLE]

and we obtain the required expression. Since $K_{j}$ are positive semi-definite, and the inverse of positive semi-definite matrices is a convex function, the entire function is convex in $\bm{s}$. $\square$∎

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

[TABLE]

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

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

[TABLE]

and we solve the mixed-integer linear programming problem:

[TABLE]

to obtain $\bm{s}_{t+1},\eta_{t+1}$. We see that $\eta_{t+1}$ is exactly the minimum value of the current approximation of $c(\bm{s})$, $c_{t}(\bm{s})$, defined below:

[TABLE]

Since $c(\bm{s})$ is convex, the piecewise linear approximation $c_{t}(\bm{s})$ is an outer approximation ($c_{t}(\bm{s})<=c(\bm{s})\;\forall\bm{s}$), so $\eta_{t}\leq c(\bm{s}_{t})\;\forall t$. As the algorithm progresses, the set of linear approximations form an increasingly better approximation of $c(\bm{s})$, and $\eta_{t}$ increases with $t$. The algorithm terminates once $\eta_{t}$ does not further increase, as it implies the linear approximation shares the same minimum value as the true function $c(\bm{s})$, which is the desired value.

Once the optimal solution $\bm{s}^{*}$ is reached, we can obtain the optimal $\bm{U}$ using the closed form solution in (6) and recover $\bm{X}$. In the next section, we discuss how this algorithm can be implemented in the context of $c(\bm{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(\bm{s})$ and its derivative $\nabla c(\bm{s})$. We first introduce the notations $\alpha_{i}(\bm{s})\in\mathbb{R}$ and $\bm{\gamma}_{i}(\bm{s})\in\mathbb{R}^{m\times 1}$.

[TABLE]

Then, the function $c(\bm{s})$ in (7) can be expressed as

[TABLE]

To calculate the derivative $\nabla c(\bm{s})$, it is easier to utilize the expression in Theorem 1 and then utilize the chain rule. After some manipulations, we obtain

[TABLE]

Therefore, we would focus on calculating $\bm{\gamma}_{i}(\bm{s})$. First, by the Matrix Inversion Lemma (Woodbury (1949)) we have

[TABLE]

where $\bm{V}\in\mathbb{R}^{m\times k}$ is the feature matrix formed by the $k$ columns of $\bm{B}$ such that $s_{j}=1$, and we have suppressed the dependency of $\bm{V}$ on $\bm{s}$ for notation ease. Note that in order to compute $\gamma_{i}(\bm{s})$ using Eq. (9) we need to invert an $m\times m$ matrix, while from Eq. (12) we need to invert a $k\times k$ matrix $\frac{\bm{I}_{k}}{\gamma}+\bm{V}^{T}\bm{W}_{i}\bm{V}$, which only requires $O(k^{3})$ calculations. Furthermore, note that calculating $\overline{\bm{a}}_{i}\overline{\bm{a}}_{i}^{T}$ only requires $|\Omega_{i}|$ multiplications where $\Omega_{i}$ is the number of known entries in row $i$ of $\bm{A}$ as we do not need to multiply on the unknown entries. Similarly, we can compute $\bm{V}^{T}\overline{\bm{a}}_{i}^{T}$ in $|\Omega_{i}|k$ multiplications, and $\bm{V}^{T}\bm{W}_{i}\bm{V}$ in $|\Omega_{i}|k^{2}$ multiplications.

Therefore, we can compute $\gamma_{i}(\bm{s})$ in floating point complexity of $O(|\Omega_{i}|k^{2}+k^{3})$. Then to calculate $\overline{\bm{a}}_{i}\bm{\gamma}_{i}(\bm{s})$ in (10) and $-\gamma\left(\bm{B}^{T}\bm{W}_{i}\bm{\gamma}_{i}(\bm{s})\right)^{2}$ (11) only requires $O(|\Omega_{i}|)$ and $O(|\Omega_{i}|p)$ calculations respectively. Thus, the total complexity of generating a full cutting plane is:

[TABLE]

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:

[TABLE]

where $\bm{U}\in\mathbb{R}^{n\times p_{1}}$ is a known matrix of $p_{1}$ features of each row, $\bm{B}\in\mathbb{R}^{m\times p_{2}}$ is a known matrix of $p_{2}$ features of each column, and $\bm{L}\in\mathbb{R}^{p_{1}\times p_{2}}$ is a sparse matrix that has $k$ nonzero entries, ensuring that $\operatorname{Rank}(\bm{X})\leq k$. We note that in Eq. (14) we restrict the support of matrix $\bm{L}$ to be $k$, rather than forcing the entries of $\bm{L}$ to be binary. This is because unlike in the one-sided case, both $\bm{U}$ and $\bm{B}$ are known, so we cannot apply the scaling transformation in (3).

We denote by $\bm{u}^{i}\in\mathbb{R}^{n\times 1}$ the $i$th column of $\bm{U}$ and $\bm{b}^{j}\in\mathbb{R}^{m\times 1}$ the $j$th column of $\bm{B}$. We introduce the matrices $\bm{W}_{i}$ as in Theorem 1. Using $\bm{X}=\bm{U}\bm{L}\bm{B}$, we can write

[TABLE]

where $D^{q,\ell}_{ij}=(\bm{u}^{q}(\bm{b}^{\ell})^{T})_{ij}$ is the $(i,j)$th entry of the matrix formed by multiplying $q$th column of $\bm{U}$ with $\ell$th column of $\bm{B}$. Then, Problem (14) becomes:

[TABLE]

As every $\bm{D}$ matrix is known, this becomes a sparse regression problem where there are $p_{1}p_{2}$ features to choose from (the $\bm{D}$ matrices), there are $|\Omega|$ samples (the $\bm{A}$ matrix), the sparsity requirement is $k$, the regression coefficients are $\bm{L}$, and we have Tikhonov regularization. Vectorizing $\bm{D}$, $\bm{L}$, and $\bm{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 $O(|\Omega|(p+k^{2})+nk^{3})$ calculations. However, in very high dimensions where $|\Omega|,n,m$ 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:

[TABLE]

where:

[TABLE]

We approximate the inner term $\alpha_{i}(\bm{s})$ by choosing $1\leq f<m$ samples from $\{1,\ldots,m\}$ without replacement, with the set denoted $F$. Then we formulate the submatrix $\bm{V}_{F}$ with such selected rows, and similarly with $\overline{\bm{a}}_{Fi}$. Then we calculate the approximation:

[TABLE]

Then we choose $1\leq g<n$ samples from $\{1,\ldots,n\}$ without replacement, with the set denoted $G$. We can then calculate an approximation of $c(\bm{s})$ using the approximated $\alpha^{F}_{i}(\bm{s})$:

[TABLE]

where the set $F$ is chosen independently for every row $i\in G$. Then the derivative of $\tilde{c}_{G}^{F}(\bm{s})$ is:

[TABLE]

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 $\tilde{c}^{F}_{G}(\bm{s})$ 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 $\bm{A}$ be a partially known matrix, $\bm{B}$ a known feature matrix, and $\bm{W}_{i}$ as defined in Theorem 1. Let $F$ be a random sample of size $f$ from the set $\{1,\ldots,m\}$, chosen without replacement. With probability at least $1-\epsilon$, we have

[TABLE]

where $M,M^{\prime}$ are absolute constants.

We see that, without assumptions on the data, the inner approximation for both the value and the derivative follows a bound with $O\left(\sqrt{\frac{(p+k)}{f}}\right)$ terms with very high probability. Furthermore, inverting the statements give that, for all $i\in\{1,\ldots,m\}$ and all $\bm{s}\in S_{k}^{p}$:

[TABLE]

So the failure probability drops off exponentially with increasing bound $\delta$, reflecting a Gaussian tail structure for $\alpha_{i}(\bm{s})$ and $\bm{B}^{T}\bm{W}_{i}\bm{\gamma}_{i}^{F}(\bm{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 $c^{F}_{G}(\bm{s})$ and $\nabla c^{F}_{G}(\bm{s})$:

Theorem 3

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

[TABLE]

where $A,B$ are absolute constants.

Similar to the inner approximations, $c^{F}_{G}(\bm{s})$ and $\nabla c^{F}_{G}(\bm{s})$ 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(\bm{s})$ is defined as the largest positive number for which the the following statement is true:

[TABLE]

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:

Proposition 1

Assume that there does not exist $\bm{s}_{1},\bm{s}_{0}\in S^{p}_{k}$ such that $c(\bm{s}_{1})=c(\bm{s}_{0})$. Then $a>0$.

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

Theorem 4

For the matrix completion problem (4), let $\bm{B}\in\mathbb{R}^{m\times p}$ be a known feature matrix, $\bm{A}\in\mathbb{R}^{n\times m}$ 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-kC\exp\left(-\frac{Da^{4}g}{(p+k)}\right)$ where $D$ is an absolute constant independent of $C,f,g,k,p$, and $a$ is the convexity parameter of the functions $\tilde{\alpha}_{i}^{s}(\bm{s})$.

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 $O\left(\exp\left(\frac{-g}{(p+k)}\right)\right)$. 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ès and Tao (2010) showed that to complete a square $N\times N$ matrix of rank $k$, we need at least $O(kN\log N)$ elements. Assume an average known rate of $\alpha=\frac{|\Omega|}{mn}$ in the original matrix $\bm{A}$, the expected number of known elements under a sampling of $f$ and $g$ is $\alpha fg$. Using $N^{2}=mn$, we need that:

[TABLE]

for some constant $c$. Theorem 4 showed that the bound on failure probability scales with $O\left(\exp\left(\frac{-g}{(p+k)}\right)\right)$, and thus we cannot have $g$ too small. Using (13), the complexity of the cutting plane with $f$ and $g$ samples are:

[TABLE]

Therefore, if we fix the expected known elements ($\alpha fg$) constant, it is more advantageous to select a smaller $g$, as $g$ scales with $k^{3}$. Thus, we set:

[TABLE]

Experimentally, we found $g_{0}=100$, $c=1$ to generate good results (the results were similar for $\frac{1}{2}\leq c\leq 2$). Therefore, by (18), the approximated cutting plane has a computational complexity of:

[TABLE]

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 $\bm{A}=\bm{UV}+\bm{E}$, where $\bm{U}\in\mathbb{R}^{n\times k}$, $\bm{V}\in\mathbb{R}^{k\times m}$, and $\bm{E}$ is an error matrix with individual elements sampled from $N(0,0.01)$. We sample the elements of $\bm{U}$ and $\bm{V}$ from a uniform distribution of $[0,1]$, and then randomly select a fraction $\mu=1-\alpha$ to be missing. We formulate the feature matrix $\bm{B}$ by combining $\bm{V}\in\mathbb{R}^{k\times m}$ with a confounding matrix $\bm{Z}\in\mathbb{R}^{(p-k)\times m}$ 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 $(m,n,p,k,\mu)$, we ran 10 tests and report the median value for every statistic.

We report the following statistics with $\bm{s}^{*}$ being the ground-truth factor vector, and $\overline{\bm{s}}$ the estimated factor vector.

• •

$n,m$ - the dimensions of $\bm{A}$.

• •

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

• •

$k$ - the true number of features.

• •

$\mu$ - The fraction of missing entries in $\bm{A}$.

• •

$T$ - the total time taken for the algorithm.

• •

MAPE - the Mean Absolute Percentage Error (MAPE) for the retrieved matrix $\hat{\bm{A}}$:

[TABLE]

where $\mathcal{S}=\Omega^{c}$ is the set of missing data in $\bm{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 $\gamma$ of OptComplete, the rank parameter of IMC and the penalization parameter $\lambda$ 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 $n,m,p,k,\mu$ 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.

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.

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 $n=10^{6}$ and $m=10^{4}$, OptComplete is over 20 times faster than IMC and over 40 times faster than SIALS. At the scale of $n=10^{6}$ and $m=10^{5}$, 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 $60$x at the scale of $n=10^{5}$ and $m=10^{4}$.

We analyze the scaling of OptComplete as a function of:

$\mu$ - 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 $|\Omega|$. 2. 2.

$n$ - The algorithm has good scalability in $n$, reflecting its $O(\sqrt{n}\log(n))$ type complexity. This allows the algorithm to support matrices with $n$ in the $10^{6}$ range. Its scaling behavior is superior to both IMC and SIALS. 3. 3.

$m$ - The algorithm scales exceptionally well in $m$. We observe that empirically the algorithm runtime seems to grow much slower than the theoretical $O(\sqrt{m}\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 $O(\sqrt{m})$ rather than the $O(m)$ as expected. We see that IMC and SIALS scales as $O(m)$. 4. 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)-O(p)$ cutting planes. Thus, as each cutting plane is $O(p)$, we expect $O(p)-O(p^{2})$ 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. 5. 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 ($k\leq 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 $\bm{A}$ where $A_{ij}$ represents the rating of individual $i$ on movie $j$, the goal is to complete the matrix $\bm{A}$ under a low-rank assumption.

The feature matrix $\bm{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:

Base - $\bm{A}_{1}$ has dimensions $3,923\times 103$. 2. 2.

Small - $\bm{A}_{2}$ has dimensions $18,227\times 323$. 3. 3.

Medium - $\bm{A}_{3}$ has dimensions $96,601\times 788$. 4. 4.

Large - $\bm{A}_{4}$ has dimensions $471,268\times 1760$. 5. 5.

Full - $\bm{A}$ has dimensions $471,268\times 14,538$.

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

For each individual matrix, we uniformly randomly withhold $20\%$ of the ratings as a test set $\mathcal{S}$, and use the remaining $80\%$ of ratings to impute a complete matrix $\hat{\bm{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:

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.

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.

Appendix A Proof of Theorem 2

We first note that since $S_{k}^{p}$ is a finite set, we only need to prove the result for a particular $\bm{s}\in S_{k}^{p}$, and it would apply for all $\bm{s}$. Therefore, we would assume $\bm{s}$ is fixed below. For simplicity, we would only demonstrate the proof for $\alpha_{i}(\bm{s})=\overline{\bm{a}}_{i}\bm{\gamma}_{i}(s)$, as the one for $(\bm{B}\bm{W}_{i}\bm{\gamma}_{i}(s))^{2}$ follows in the same exact fashion. Furthermore, since we are only focusing on one particular $i\in\{1,\ldots,n\}$, we would drop all $i$ subscripts below for ease of notation. The quantities of interest are therefore:

[TABLE]

First, let $\overline{\bm{V}}=\bm{W}\bm{V}$, and let us consider a reduced QR factorization of $\overline{\bm{V}}=\bm{Q}\bm{R}$ where $\bm{Q}\in\mathbb{R}^{m\times k}$ has orthogonal columns such that $\bm{Q}^{T}\bm{Q}=m\cdot\bm{I}_{k}$, and $\bm{R}_{i}\in\mathbb{R}^{k\times k}$. Note such definition implies $\|\bm{R}\|=O(1)$. Then, we would rewrite the terms as follows:

[TABLE]

We note that

[TABLE]

Therefore, if we treat $\overline{a}_{1}^{2},\ldots\overline{a}_{m}^{2}$ as a finite population, then $\overline{\bm{a}}_{F}\overline{\bm{a}}_{F}^{T}$ is a random sample of $f$ points drawn without replacement from that set, and similarly for $\overline{\bm{a}}_{F}\bm{V}_{F}$. Therefore, we can then utilize Hoeffding’s inequality to bound the deviation of these terms, as reproduced below:

Proposition 2* *(Hoeffding’s Inequality)

Let $\mathcal{X}=(x_{1},\ldots,x_{n})$ be a finite population of $N$ points and $X_{1},\ldots,X_{n}$ be a random sample drawn without replacement from $\bm{X}$. Let

[TABLE]

Then, for all $\epsilon>0$, we have

[TABLE]

For a proof, see for example Boucheron et al. (2013). Then, applying Proposition 2 to $\overline{\bm{a}}_{f}\overline{\bm{a}}_{f}^{T}$, $\overline{\bm{a}}_{f}\bm{V}_{f}$, and inverting the inequality, we have

[TABLE]

where $A,B$ are constants independent of $k,f,m,\epsilon$.

Now we would show that $\frac{\bm{R}^{T}\bm{Q}^{T}\bm{Q}\bm{R}}{m}$ is close to $\frac{\bm{R}^{T}\bm{Q}_{F}^{T}\bm{Q}_{F}\bm{R}}{f}$:

Lemma 1

[TABLE]

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):

Lemma 2

Let $\mathcal{X}\in\mathbb{R}^{k\times k}$ be a finite set of positive-semidefinite matrices, and suppose that

[TABLE]

where $\lambda_{\min}/\lambda_{\max}$is the minimum/maximum eigenvalue function. Sample $\{\bm{X}_{1},\ldots,\bm{X}_{\ell}\}$ uniformly at random without replacement. Compute:

[TABLE]

Then,

[TABLE]

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

Proof.

Proof of Lemma 1: First, we note that

[TABLE]

where $\bm{q}_{i}^{T}\bm{q}_{i}\in\mathbb{R}^{k\times k}$ rank-one positive semi-definite matrices. Therefore, we can take $\bm{Q}^{T}_{F}\bm{Q}_{F}$ as a random sample of size $f$ from the set $\mathcal{X}=\{\bm{q}_{i}^{T}\bm{q}_{i}\}_{i=1,\ldots,m}$, which satisfies the conditions in Lemma 2 with $D=O(k)$. Furthermore, with $\mathcal{X}$, we observe that we have $\mathbb{E}\bm{X}_{1}=\frac{\bm{Q}^{T}\bm{Q}}{m}=\bm{I}_{k}$, so we have

[TABLE]

Therefore, we apply Lemma 2 to $\bm{Q}_{F}^{T}\bm{Q}_{F}$ and obtain

[TABLE]

where we set $D=\frac{kD^{\prime}}{4}$ with $D^{\prime}=O(1)$. Some rearrangement gives:

[TABLE]

Now since $\frac{\bm{Q}^{T}\bm{Q}}{m}=\bm{I}_{k}$, we have

[TABLE]

Combining equation (A6) and (A5) gives

[TABLE]

Then, we have

[TABLE]

Taking $C=D^{\prime}\|\bm{R}\|^{4}\log(2)$ gives the required result. (Note $\|\bm{R}\|=O(1)$ as we setup the QR decomposition to have $\bm{Q}^{T}\bm{Q}=O(m)$). $\square$∎

With equations (A2), (A3), and (A4), we are now ready to bound $\alpha(\bm{s})$ and $\alpha^{F}(\bm{s})$. We first introduce another lemma from matrix perturbation theory (for proof, see e.g. Stewart (1990)).

Lemma 3

Let $\bm{A},\bm{B}$ be invertible matrices and let $\bm{B}=\bm{A}+\bm{\Delta}$. Then, we have

[TABLE]

Then, we have

[TABLE]

For some $O(1)$ constant $M$ with probability at least $1-4\epsilon$. Therefore, taking $\epsilon^{\prime}=\frac{\epsilon}{4}$ and absorbing the extra constant into $M$ gives the requires result.

Appendix B Proof of Theorem 3

We first prove a Hoeffding-type bound as follows

Proposition 3

Let $X_{1},\ldots,X_{n}$ be independent (but not necessarily identically distributed) random variables which satisfy

[TABLE]

Then, we have

[TABLE]

We note that such statement differs from Hoeffding’s inequality as we do not require $\mathbb{E}[X_{i}]=a_{i}$. The proof is as follows.

Proof.

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

Lemma 4

Let $Z$ be a random variable on $\mathbb{R}$. Then for all $t>0$, we have

[TABLE]

where $M_{Z}(s)$ is the moment generating function of $Z$.

Let $Y_{i}=X_{i}-a_{i}$. Then, we have

[TABLE]

We aim to bound $\mathbb{E}\left[e^{sY_{i}}\right]$, subject to $\mathbb{P}(|Y_{i}|\geq t)\leq\exp\left(\frac{-t^{2}}{\sigma_{i}^{2}}\right)$. Since $e^{st}$ is an increasing function of $t$, $\mathbb{E}\left[e^{sY_{i}}\right]$ is maximized when $\mathbb{P}(Y_{i}\geq t)=\exp\left(\frac{-t^{2}}{\sigma_{i}^{2}}\right)$, which results in a probability distribution function of $f_{Y_{i}}(y)=2\frac{t}{\sigma_{i}^{2}}\exp\left(\frac{-t^{2}}{\sigma_{i}^{2}}\right)$ for $Y_{i}$. Then, we have

[TABLE]

We then substitute this result into (A14) and obtain

[TABLE]

Therefore,

[TABLE]

By applying the previous derivation to $-Y_{1},\ldots,-Y_{n}$, we obtain

[TABLE]

Combining the two results completes the proof. $\square$∎

Using Proposition 3, we bound the difference between $\tilde{c}^{F}_{G}(\bm{s}_{t})$ and $c(\bm{s}_{t})$. We have

[TABLE]

with probability $1-2\epsilon$. Therefore, taking $\epsilon^{\prime}=2\epsilon$ and $A=M^{\prime\prime}\log(2)$ we have

[TABLE]

Through a similar derivation, we have

[TABLE]

Appendix 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 $\bm{T}_{ij}=\frac{\bm{I}}{\gamma}+\bm{W}_{i}\bm{K}_{j}\bm{W}_{i}$. Then our cost function can be rewritten as:

[TABLE]

Then the second derivative can be easily calculated as:

[TABLE]

Thus, for any vector $\bm{t}\in\mathbb{R}^{p}$, we have:

[TABLE]

Now define $\bm{v}_{i}=\left(\sum_{l=1}^{p}t_{l}\bm{T}_{il}\right)\left(\sum_{j=1}^{p}s_{j}\bm{T}_{ij}\right)^{-1}\bar{\bm{a}}_{i}^{T}$, then our expression becomes:

[TABLE]

Now since $a=0$, the Hessian cannot be positive definite for all $\bm{s}$ (as if it was, then $a>0$). Then there exist $\bm{s}^{0},\bm{s}^{1}\in S^{p}_{k}$ such that for $\bm{t}=\bm{s}^{0}-\bm{s}^{1}$, we have:

[TABLE]

and $\bm{t}$ Note that for any $\gamma>0$, $\bm{T}_{ij}$ is positive definite, so $\sum_{j=1}^{p}s_{j}\bm{T}_{ij}$ is positive definite for any $\bm{s}\in S^{p}_{k}$ and all $i$. Therefore if $a=0$, we must have $\bm{v}_{i}=0$ for all $i$, which means that:

[TABLE]

for all $i$. Now note that we have:

[TABLE]

We have that:

[TABLE]

And similarly, we have \bm{t}^{T}\frac{\partial c(\bm{s})}{\partial s_{k}\partial s_{l}}\bigg{\rvert}_{\bm{s}=\bm{s}^{0}}\bm{t}=0 and all higher derivatives being 0.

Thus, we have:

[TABLE]

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$.

Appendix 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 $\bm{s}^{*}$ be an optimal solution for Problem (4). Let $\bm{s}_{t}$ be an optimal solution at the $t$-th iteration of OptComplete, $t\in[C]$. The cutting plane constraint for OptComplete at the point of an optimal solution $\bm{s}^{*}$ is

[TABLE]

If $c(\bm{s}^{*})<\eta_{t}$, then $\bm{s}^{*}$ will be cut off, and OptComplete will not find $\bm{s}^{*}$. Applying the definition of the convexity parameter (16) and letting $\|\bm{s}^{*}-\bm{s}_{t}\|=\theta_{t}$ (noting that $\theta_{t}\geq 1$) we obtain

[TABLE]

Therefore, if

[TABLE]

or equivalently if

[TABLE]

then OptComplete will not find $\bm{s}^{*}$. Therefore, for OptComplete to succeed, $\zeta_{t}$ should satisfy $\zeta_{t}\leq\theta_{t}^{2}a^{2}/2$ for all $t\in[C].$

Let $S=$ the event that OptComplete succeeds in finding $\bm{s}^{*}$. Then,

[TABLE]

Since at each step of OptComplete we randomly sample $r$ new rows and $s$ new columns, the events $\zeta_{t}\leq\frac{\theta_{t}^{2}a^{2}}{2}$ are independent for different $t\in[C]$, and hence

[TABLE]

Therefore, we focus on calculating $\mathbb{P}(\zeta_{t}\geq\theta_{t}^{2}a^{2}/2)$. Since

[TABLE]

using Theorem 3, we can thus provide the following bound for deviation of $\zeta_{t}$, for some constant $C$:

[TABLE]

where $\theta_{t}=\|\bm{s}^{*}-\bm{s}_{t}\|$. Then we can invert this to calculate $\mathbb{P}(\zeta_{t}\geq\theta_{t}^{2}a^{2}/2)$

[TABLE]

taking $D=\frac{1}{4C}$, and noting that $\theta_{t}\geq 1$. Then, we substitute (A25) into (A23) to obtain

[TABLE]

completing the proof.

Appendix E List of Features Used in the Netflix 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 multicollinearity)

• •

Number of Original Music Score

• •

Number of Male Actors

• •

Number of Female Factors

• •

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

• •

Constant