Convex relaxations of structured matrix factorizations

Introduction

Structured matrix factorization has many applications in various areas of science engineering, i.e., clustering and principal component analysis , source separation , signal processing , machine learning , and all domains where reduced representations are desired.

Without any structure, traditional principal component analysis may be solved exactly in polynomial time through a singular value decomposition. However, adding additional structure on the components $U$ and $V$ of the factorization of $X=UV^{\top}$ (e.g., sparsity, non-negativity or discreteness) is most often done through alternating minimization (with respect to $U$ and $V$). While all steps are usually done through convex optimization, the problem is not jointly convex, and there are typically multiple local minima, and algorithms usually come with no convergence guarantees. For example, in presence of positivity constraints, the problem of non-negative matrix factorization (NMF) may not be solved in polynomial-time in general, and most algorithms for NMF (e.g., ) perform a form a block coordinate descent with no guarantees (see hardness results and particular situations of actual solvability in ). In this paper, we follow a convex relaxation approach.

We impose some structure on each column of $U$ and $V$ and consider a general convex framework which amounts to computing a certain gauge function (such as a norm) at $X$, from which the decomposition may be obtained (for example, the nuclear norm leads to the usual singular value decomposition). This convex framework corresponds to removing any rank constraint on the decomposition and has appeared under various forms in the literature, as summing norms , decomposition norms or (a special case of) atomic norms . This is presented in details (equivalent representations, rotation-invariant cases, weighted nuclear norm formulations) in Section 4. An interesting aspect is that the gauge functions we consider are polar to generalizations of matrix norms, which are commonly used in many areas of applied mathematics, in particular in robust optimization and control .

The statistical and recovery properties of these norms and their relaxations have been studied in several contexts ; in this paper, we focus on optimization aspects. The convex framework we introduce in Section 4 only lead to polynomial-time algorithms in few situations (e.g., the nuclear norm based on the singular value decomposition). In Section 5, we consider computable additional relaxations based on semi-definite programming. These may be used to compute the related gauge functions as well as their polar, with constant-factor approximation guarantees in some cases (see Section 5.1). The first setting where one can get dimension-independent guarantees has already been studied by and corresponds to gauge functions that have variationl diagonal representations. We also consider a more general setting with dimension-dependent bounds.

A key practical problem is to obtain not only a lower-bound on the value of the gauge function, but also an explicit decomposition which preserves the approximation guarantees. We present in Section 6 iterative conditional gradient algorithms and their analysis, which extend existing results in several ways: (a) we obtain convergence guarantees even when the polar gauge function may be approximately computed with a multiplicative approximatio ratio—earlier work considers only additive approximations, (b) following , they may be applied to penalized versions of the problem, i.e., to solve a generalized basis pursuit problem, (c) under some additional assumptions, they may find approximate decompositions of $X$ that converge linearly.

Finally, our framework relies on variational representations of gauge functions and in particular norms as maxima and minima of quadratic functions. These representations have been already used in several contexts (machine learning, signal processing, optimization, see, e.g., and references therein). In this paper, we provide in Section 3 a thorough analysis of these decompositions (minimal and maximal representations, duality between lower and upper bounds, sufficient and necessary conditions for diagonal or rotation-invariant representations).

Review of gauge function theory

In this section, we present relevant concepts and results from convex analysis. These tools are needed because the type of structure we want to impose go beyond what can be characterized by norms (such as positivity). See for more details on gauge functions and their properties.

Polar sets and functions.

Given any set ${\mathcal{C}}$ (not necessarily convex), the polar of ${\mathcal{C}}$ is the set ${\mathcal{C}}^{\circ}$ defined as

It is always closed and convex. Moreover, the polar of ${\mathcal{C}}$ is equal to the polar of the closure of ${\rm hull}({\mathcal{C}}\cup\{0\})$.

When ${\mathcal{C}}$ is the unit ball of a norm $\Omega$, ${\mathcal{C}}^{\circ}$ is the unit ball of the dual norm which we denote $\Omega^{\circ}$ (instead of the usual definition $\Omega^{\ast}$, because the Fenchel conjugate of $\Omega$ is not the dual norm, but the indicator function of the dual unit ball).

If ${\mathcal{C}}$ is a closed convex set containing the origin, then ${\mathcal{C}}^{\circ\circ}={\mathcal{C}}$—more generally, for any set ${\mathcal{C}}$, ${\mathcal{C}}^{\circ\circ}$ is the closure of ${\rm hull}({\mathcal{C}}\cup\{0\})$. The polarity is a one-to-one mapping from closed convex sets containing the origin to themselves. In this paper, we will also consider gauge functions associated with closed potentially non convex sets ${\mathcal{C}}$, in which case, we mean the gauge function associated to ${\mathcal{C}}^{\circ\circ}={\rm hull}({\mathcal{C}}\cup\{0\})$, i.e., $\gamma_{\mathcal{C}}=\gamma_{{\mathcal{C}}^{\circ\circ}}$.

Operations on gauge functions.

Links with convex hulls.

Given a compact set ${\mathcal{P}}$ and its compact convex hull ${\mathcal{C}}$ (for example, ${\mathcal{P}}$ might be the set of extreme points of ${\mathcal{C}}$), we have ${\mathcal{P}}^{\circ}={\mathcal{C}}^{\circ},$ since maxima of linear functions on ${\mathcal{C}}$ or ${\mathcal{P}}$ are equal. An alternative definition of $\gamma_{\mathcal{C}}$ is then

Moreover, in the definition above, by Caratheodory’s theorem for cones, we may restrict the cardinality of $I$ to be less than or equal to $d$.

Representations of gauge functions through quadratic functions

We first consider closed convex sets ${\mathcal{K}}\subset{\mathcal{S}}_{d}^{+}$ of symmetric positive definite matrices such that

When ${\mathcal{U}}$ is symmetric, the variational formulation in Eq. (1) leads to a representation of $\gamma_{\mathcal{U}}(x)^{2}$ as a convex function $I_{\mathcal{K}}^{\ast}(xx^{\top})$ of $xx^{\top}$. Note that in general, ${\mathcal{K}}$ is not unique. We now show that there always exists a set ${\mathcal{K}}$ satisfying Eq. (1), and provide a description of the largest such set.

Let ${\mathcal{P}}_{\mathcal{U}}=\{xx^{\top},\ x\in{\mathcal{U}}\}$. Then ${\mathcal{P}}_{\mathcal{U}}^{\circ}\cap{\mathcal{S}}_{d}^{+}$ is the larget closed convex set ${\mathcal{K}}\subset{\mathcal{S}}_{d}^{+}$ of positive semidefinite matrices such that Eq. (1) is satisfied.

We now have $\sup_{M\in{\mathcal{P}}_{\mathcal{U}}^{\circ}\cap{\mathcal{S}}_{d}^{+}}\sqrt{x^{\top}Mx}\leqslant\sup_{M\in{\mathcal{P}}_{\mathcal{U}}^{\circ}}\sqrt{x^{\top}Mx}=\min\{\gamma_{\mathcal{U}}(x),\gamma_{\mathcal{U}}(-x)\}$. Since this is a convex function of $x$, it must be less than its convex envelope, which is $\gamma_{{\mathcal{U}}\cup(-{\mathcal{U}})}(x)$ (since they have the same Fenchel conjugates).

Moreover, if $v\in{\mathcal{U}}^{\circ}\cap(-{\mathcal{U}})^{\circ}$, then $vv^{\top}\in{\mathcal{P}}_{\mathcal{U}}^{\circ}$. This implies that $\sup_{M\in{\mathcal{P}}_{\mathcal{U}}^{\circ}}x^{\top}Mx\geqslant\sup_{v\in{\mathcal{U}}^{\circ}\cap(-{\mathcal{U}})^{\circ}}(v^{\top}x)^{2}=I^{\ast}_{{\mathcal{U}}^{\circ}\cap(-{\mathcal{U}})^{\circ}}(x)^{2}=\gamma_{{\mathcal{U}}^{\circ}\cap(-{\mathcal{U}}^{\circ})}^{\circ}(x)^{2}=\gamma_{{\mathcal{U}}\cup(-{\mathcal{U}})}(x)^{2}$. Thus Eq. (1) is indeed satisfied by ${\mathcal{P}}_{\mathcal{U}}^{\circ}\cap{\mathcal{S}}_{d}^{+}$. Finally, if ${\mathcal{K}}$ satisfies Eq. (1), then we must have ${\mathcal{K}}\subset{\mathcal{P}}_{\mathcal{U}}^{\circ}$ by definition of polar sets, hence ${\mathcal{P}}_{\mathcal{U}}^{\circ}\cap{\mathcal{S}}_{d}^{+}$ is the largest.

Note that the set ${\mathcal{P}}_{\mathcal{U}}^{\circ}\cap{\mathcal{S}}_{d}^{+}$ is equal to $\{M\in{\mathcal{S}}_{d}^{+},\ \forall u\in{\mathcal{U}},u^{\top}Mu\leqslant 1\}$—this representation was already considered in for norms. In certain situations, the largest possible set is desirable (for example when deriving convex relaxations). In other situations (for example when using these representations for optimization), smallest sets are desirable. However, such a notion is not possible. Indeed, for $\gamma_{\mathcal{U}}=\|\cdot\|_{2}$, the sets ${\mathcal{K}}=\{I\}$ and ${\mathcal{K}}=\{M\in{\mathcal{S}}_{d}^{+},\ \|M\|_{F}\leqslant 1\}$, are two possible sets, and thus there is no single smallest set. One possible small set is the convex hull of the maximal elements of ${{\mathcal{P}}}_{\mathcal{U}}^{\circ}\cap{\mathcal{S}}_{d}^{+}$ (for the positive semi-definite order).

In the proof of Prop. 1, we have introduced the gauge function $\gamma_{{\mathcal{P}}_{\mathcal{U}}}$. We now provide a representation of $\gamma_{{\mathcal{P}}_{\mathcal{U}}}$ related to factorizations of positive semidefinite matrices (note that in the following proposition, we only assume that ${\mathcal{U}}$ is closed and contains 0 in its hull).

Let ${\mathcal{P}}_{\mathcal{U}}=\{xx^{\top},\ x\in{\mathcal{U}}\}\subset{\mathcal{S}}_{d}^{+}$. We have, for all positive semidefinite matrix $M\in{\mathcal{S}}_{d}^{+}$:

Moreover, we may choose $r\leqslant d(d+1)/2$.

Proof This is a direct application of the representation of gauge functions and and the property $\gamma_{{\mathcal{P}}_{\mathcal{U}}}(xx^{\top})=\min\{\gamma_{\mathcal{U}}(x),\gamma_{\mathcal{U}}(-x)\}^{2}$, that was shown in the proof of Prop. 1.

The last proposition provides a structured decomposition framework for positive semi-definite matrices, that will be considered for rectangular matrices in Section 4. Obtaining explicitly the decomposition $M=UU^{\top}$ from $M$ may be done with the iterative algorithms presented in Section 6. Note that by considering a representation of $U$ as $U=M^{1/2}S$ where $SS^{\top}=I$, computing $\gamma_{{\mathcal{P}}_{\mathcal{U}}}$ may be seen as a factorization problem with two factors (which can then be used in alternating minimization procedures).

2 Minima of quadratic functions

We now consider closed convex sets ${\mathcal{L}}\subset{\mathcal{S}}_{d}^{+}$ such that

Here, we define $x^{\top}M^{-1}x$ as $x^{\top}M^{-1}x=\inf t$ such that \Big{(}\!\begin{array}[]{cc}M&x\\ x^{\top}&t\end{array}\!\Big{)}\succcurlyeq 0\Leftrightarrow tM\succcurlyeq xx^{\top}. This implies that the value may be finite even when $M$ is not invertible.

When ${\mathcal{U}}$ is symmetric, the variational formulation in Eq. (2) leads to a representation of $\gamma_{\mathcal{U}}(x)^{2}$ as a concave $\inf_{M\in{\mathcal{L}}}\mathop{\rm tr}M^{-1}xx^{\top}$ of $xx^{\top}$. This is to be contrasted with the fact that it is also a convex function of $xx^{\top}$ because of the representation discussed in Section 3.1. The two properties are in fact related through a duality argument:

Let ${\mathcal{L}}\subset{\mathcal{S}}_{d}^{+}$ be a closed convex set. Then the following two properties are equivalent:

This implies that the largest set ${\mathcal{L}}$ such that (a) is valid is ${\mathcal{P}}^{\circ}_{{\mathcal{U}}^{\circ}}\cap{\mathcal{S}}_{d}^{+}$ defined in Prop. 1.

Similarly, if we have a concave upper-bound of $\max\{\gamma_{\mathcal{U}}^{\circ}(y),\gamma_{\mathcal{U}}^{\circ}(-y)\}^{2}$ based on ${\mathcal{K}}$, we also have a convex upper-bound based on ${\mathcal{K}}^{\circ}\cap{\mathcal{S}}_{d}^{+}$. Note that it is not true in general that exact representations of one kind transfer to exact representation of the other kind.

3 Examples

4 Diagonal representations

We now consider cases where the set of matrices ${\mathcal{K}}$ and ${\mathcal{L}}$ are diagonal, i.e., all principal axes of ellipsoids are aligned with the canonical basis. For simplicity, we consider only sets ${\mathcal{U}}$ which are compact, have zero in their interior, and are invariant by sign flips of any components. The corresponding gauge functions are then absolute norms, which are functions of the absolute values of each component . The following proposition provides several characterizations (see related work in ).

Proof It is straightforward to see that (a) implies both (b) and (c).

5 Matrix norms invariant by rotation

$W\mapsto\Omega(W)^{2}$ is a concave function of $WW^{\top}$.

$V\mapsto\Omega^{\circ}(V)^{2}$ is a concave function of $VV^{\top}$.

Gauge functions and structured matrix factorizations

This function is exactly the gauge function of the set $\{uv^{\top},\ u\in{\mathcal{U}},v\in V\}$. It is therefore convex, non-negative and positively homogeneous. Moreover, the minimization problem defining it is attained for $r\leqslant nd$. There are several equivalent formulations which we are going to use, which relies on the fact that $u_{m}v_{m}^{\top}=(u_{m}\lambda_{m})(v_{m}\lambda_{m}^{-1})^{\top}$ for any $\lambda_{m}>0$ (in the following expressions, we always minimize with respect to the rank and thus omit the notation $\inf_{r\geqslant 0}$).

Moreover, $\Theta$ is a norm as soon as ${\mathcal{U}}$ and ${\mathcal{V}}$ are norm balls—in this case, these norms were studied in several settings . The next proposition shows that the polar of $\Theta$ has a simple form (the proof is straightforward from the gauge function interpretation). It is a matrix norm when ${\mathcal{U}}$ and ${\mathcal{V}}$ are norm balls; moreover, in all cases, since ${\mathcal{U}}$ and ${\mathcal{V}}$ are assumed compact, $\Theta^{\circ}$ has full domain.

2 Examples

$\gamma_{\mathcal{V}}=\|\cdot\|_{2}$ and $\gamma_{\mathcal{U}}=\sqrt{\nu\|\cdot\|_{2}^{2}+(1-\nu)\|\cdot\|_{1}^{2}}$: this is a convex relaxation of sparse coding , where a decomposition with a sparse factor $U$ and low rank (small number of columns for $U$ and $V$) is looked for.

$\gamma_{\mathcal{V}}=\|\cdot\|_{2}$ and $\gamma_{\mathcal{U}}=\|\cdot\|_{4/3}$: we have $\gamma_{\mathcal{U}}^{\circ}=\|\cdot\|_{4}$, and the dual norm is such that $\Theta^{\circ}(Y)^{4}=\max_{\|v\|_{2}\leqslant 1}\sum_{i}(Yv)_{i}^{4}$, and corresponds to maximum kurtosis projections .

3 Special case: γ𝒱=∥⋅∥2\gamma_{\mathcal{V}}=\|\cdot\|_{2}

Computing the polar of $\Theta$ then corresponds to a quadratic maximization problem:

We also have a variational quadratic representation corresponding to a rotation-invariant gauge function (Section 3.5); following Section 3.1, we denote by ${\mathcal{P}}_{\mathcal{U}}$ the set $\{uu^{\top},\ u\in{\mathcal{U}}\}$. There is an explicit representation of $\Theta$ and $\Theta^{\circ}$ in terms of ${\mathcal{P}}_{\mathcal{U}}^{\circ\circ}$ (the convex hull of ${\mathcal{P}}_{\mathcal{U}}\cup\{0\}$).

Note that $\Theta$ may not have full domain.

From the earlier representation, we get \displaystyle\Theta(X)=\inf_{M\succcurlyeq 0}\big{\{}\frac{1}{2}\gamma_{{\mathcal{P}}_{\mathcal{U}}}(M)+\frac{1}{2}\mathop{\rm tr}X^{\top}M^{-1}X\big{\}}, by simply optimizing over the half line generated by $M$. We also have a representation involving all decompositions of $X$ as $UV^{\top}$:

which is straightforward from the definitions and interpretations of $\Theta$ and $\gamma_{{\mathcal{P}}_{\mathcal{U}}}$ as a gauge function (Prop. 2). In the next section, we show how to compute $U$ and $V$ from the solution of a certain convex problem—note however that this does not lead to the optimal $U$ and $V$ in Eq. (3). See more details in Section 6.

4 General case

We now consider all possible cases, beyond $\gamma_{\mathcal{V}}=\|\cdot\|_{2}$ (where $\gamma_{{\mathcal{P}}_{\mathcal{V}}}(M)=\mathop{\rm tr}M$). We only assume that ${\mathcal{U}}$ and ${\mathcal{V}}$ are compact and contain zero in their hulls. It is tempting to consider an extension of Eq. (5), by considering

Let $\widetilde{\Theta}$ be defined in Eq. (6). Then, $\widetilde{\Theta}$ is a gauge function (convex, positively homogeneous, and non-negative), and its polar has the expression:

because $\sup_{UU^{\top}=I,\ VV^{\top}=I}\mathop{\rm tr}Z^{\top}UV^{\top}=\|Z\|_{\ast}$. Finally, we have for some $u_{m},v_{m}$, \Theta(X)=\frac{1}{2}\sum_{m}\big{\{}\gamma_{\mathcal{U}}(u_{m})^{2}+\gamma_{\mathcal{V}}(v_{m})^{2}\big{\}}. Thus $\Theta(X)\geqslant\widetilde{\Theta}(X)$ because of Prop. 2, which implies the reverse inequality for polars.

We may now give several equivalent expressions for $\widetilde{\Theta}$ and $\widetilde{\Theta}^{\circ}$ which are obtained using usual convex duality.

Let $\widetilde{\Theta}$ be defined in Eq. (6). We have:

Proof The four expressions are equivalent by taking polars or using Lagrangian duality for the semi-definite cone. The only element left to prove is that any of these four statements are correct. We have:

which leads to desired result from Prop. 9. Note that we may obtain directly a relationship between $\widetilde{\Theta}^{\circ}$ and $\Theta^{\circ}$ as follows:

that is a variational formulation of $\widetilde{\Theta}$ through nuclear norms (note again that ${\mathcal{P}}_{\mathcal{U}}^{\circ}$ and ${\mathcal{P}}_{\mathcal{V}}^{\circ}$ may not be compact and $\widetilde{\Theta}(X)$ potentially infinite). This is to be contrasted with the following representation:

which is a variational upper-bound as a weighted Frobenius norm, which is not computable as a convex program. Instead, alternate minimization could be used, each of the steps being doable as a convex program.

Given the primal/dual solutions $(M,N)$, $(Q,S)$ of the dual convex optimization problems defining $\widetilde{\Theta}$ in Prop. 10, we have the following optimality conditions: Assume $Q=AA^{\top}$ with $A$ full rank (i.e., $A^{\top}A$ invertible), and $S=BB^{\top}$ with $B$ full rank. Then an optimal $Z$ is $Z=AGH^{\top}B^{\top}$, with $A^{\top}XB=G\mathop{\rm Diag}(s)H^{\top}$ is the singular value decomposition of $A^{\top}XB$. We have

which implies (since $A$ and $B$ have full rank):

This leads to $A^{\top}MA=A^{\top}XBHG^{\top}=G\mathop{\rm Diag}(s)G^{\top}$. We then take as candidates $U=(AA^{\top})^{\dagger}AG\mathop{\rm Diag}(s)^{1/2}$ and $V=(BB^{\top})^{\dagger}BH\mathop{\rm Diag}(s)^{1/2}$, leading to $UU^{\top}=M$ and $VV^{\top}=N$. Moreover, $X=UV^{\top}$ and $\widetilde{\Theta}(X)=\|A^{\top}XB\|_{\ast}=\mathop{\rm tr}\mathop{\rm Diag}(s)$. Finally, by optimality of $Q$, we have $\gamma_{{\mathcal{P}}_{\mathcal{U}}}(UU^{\top})=\mathop{\rm tr}UU^{\top}Q=\mathop{\rm tr}A^{\top}MA=\mathop{\rm tr}\mathop{\rm Diag}(s)$ and $\gamma_{{\mathcal{P}}_{\mathcal{V}}}(VV^{\top})=\mathop{\rm tr}VV^{\top}S=\mathop{\rm tr}B^{\top}NB=\mathop{\rm tr}\mathop{\rm Diag}(s)$, and thus $(U,V)$ is an optimal decomposition for Eq. (6).

The duality properties presented in Prop. 10 are valid for all convex sets of positive semi-definite matrices that contain zero (and in particular for the ones in Section 5). This will be used in Section 5 to get computable approximations (with guarantees).

5 Non-convex approaches

Finding approximations of the gauge function $\Theta$ or its polar $\Theta^{\circ}$ has been tackled thoroughly through non-convex approaches, which may only find stationary points of the associated optimization problems. In this paper, we briefly review two approaches based on alternating minimization.

In Section 6, obtaining approximate maximizers $(u,v)$ will be key for the generalized conditional gradient algorithm and simplicial methods. In our experiments in Section 7, the power method leads to better values once initialized from results of the convex relaxations presented in Section 5.

Alternating optimization.

We consider $r\leqslant np$ elements and minimize either the sum $\sum_{m=1}^{r}\gamma_{\mathcal{U}}(u_{m})\gamma_{\mathcal{V}}(v_{m})$ or \sum_{m=1}^{r}\big{\{}\frac{1}{2}\gamma_{\mathcal{U}}(u_{m})^{2}+\frac{1}{2}\gamma_{\mathcal{V}}(v_{m})^{2}\big{\}}, subject to $X=\sum_{m=1}^{r}u_{m}v_{m}^{\top}$. This may be done by alternating minimization, with no guarantees, beyond decreasing the value of the upper-bound on the norm.

Non-convex optimization without local minima.

Nevertheless, in practice, alternating optimization techniques, although they come with no guarantees, still perform well when they can be applied (see Section 7).

Computable convex relaxations of decomposition norms

Before describing computable relaxations, we first refine the computability requirement which is weaker than having a good representation of ${\mathcal{P}}_{\mathcal{U}}^{\circ\circ}$ or its polar ${\mathcal{P}}_{\mathcal{U}}^{\circ}$. We start with a simple lemma:

Let $\mathcal{A},\mathcal{B}\subset{\mathcal{S}}_{n}^{+}$. The following statements are equivalent:

(a) $\forall M\in{\mathcal{S}}_{n}^{+}$, $\sup_{N\in\mathcal{A}}\mathop{\rm tr}MN\leqslant\sup_{N\in\mathcal{B}}\mathop{\rm tr}MN$,

(b) $\mathcal{A}^{\circ\circ}-{\mathcal{S}}_{n}^{+}\subset\mathcal{B}^{\circ\circ}-{\mathcal{S}}_{n}^{+}$,

(c) $\mathcal{A}^{\circ}\cap{\mathcal{S}}_{n}^{+}\supset\mathcal{B}^{\circ}\cap{\mathcal{S}}_{n}^{+}$.

Moreover, if $\mathcal{A}$ is composed of rank-one matrices, then they are also equivalent to

Proof (a) is equivalent to (c) by simple properties of gauge functions. (b) is equivalent to (c) because of polar calculus and $({\mathcal{S}}_{n}^{+})^{\circ}=-{\mathcal{S}}_{n}^{+}$, which imply $(\mathcal{A}^{\circ}\cap{\mathcal{S}}_{n}^{+})^{\circ}={\rm hull}(\mathcal{A}^{\circ\circ}\cup(-{\mathcal{S}}_{n}^{+}))=\mathcal{A}^{\circ\circ}-{\mathcal{S}}_{n}^{+}$. (a) trivially implies (d). We only need to show (d) implies (a). Assume (d) is true and let $M=YY^{\top}\in{\mathcal{S}}_{n}^{+}$, we have

Given the representation of $\widetilde{\Theta}$ in Prop. 10, a key consequence of the previous lemma is that in the representation of $\widetilde{\Theta}(X)$ and $\widetilde{\Theta}^{\circ}(Y)$ in Section 4, we may replace ${\mathcal{P}}_{{\mathcal{U}}}^{\circ\circ}$ by any set ${\mathcal{C}}_{\mathcal{U}}$ such that

For all other cases, we need computable convex approximations ${\mathcal{C}}_{\mathcal{U}}$ of ${\mathcal{P}}_{\mathcal{U}}^{\circ\circ}$ such that for

(if we impose the property above for all symmetric matrices $M$, it is equivalent to ${\mathcal{P}}_{\mathcal{U}}^{\circ\circ}\subset{\mathcal{C}}_{\mathcal{U}}$). Given Lemma 1, this is equivalent to

that is, we need an upper-bound on $\max\{\gamma_{\mathcal{U}}^{\circ}(y),\gamma_{\mathcal{U}}^{\circ}(-y)\}^{2}$, which is convex in $yy^{\top}$. Note that given two sets ${\mathcal{C}}_{\mathcal{U}}$ that satisfy Eq. (32), their intersection also does and the relaxation is then always tighter. We can therefore use several properties of the set ${\mathcal{U}}$, that may then be combined together:

Linear inequalities: any matrix $M$ such that $\max_{\gamma_{\mathcal{U}}(u)\leqslant 1}u^{\top}Mu$ may be computed or efficiently upper-bounded by $h$ adds another constraint of the form $\mathop{\rm tr}MU\leqslant h$. Typically, $M=I$ is a good candidate, as this defines the equivalence between the gauge function $\gamma_{\mathcal{U}}$ and $\|\cdot\|_{2}$.

In Table 1, we present relaxations for the examples we consider in this paper. For all relaxations we consider (and also with positivity constraints), we have: ${\mathcal{C}}_{\mathcal{U}}\cap\{M,\ {\rm rank}(M)=1\}={\mathcal{P}}_{\mathcal{U}}$. This implies that for all $u\in{\rm dom}(\gamma_{\mathcal{U}})$, $\gamma_{{\mathcal{P}}_{\mathcal{U}}}(uu^{\top})=\gamma_{{\mathcal{C}}_{\mathcal{U}}}(uu^{\top})=\min\{\gamma_{\mathcal{U}}(u),\gamma_{\mathcal{U}}(-u)\}^{2}$. Moreover, for all relaxations in Table 1 (without positivity constraints), we have equality in Eq. (32), a property that will be useful when deriving performance guarantees in the next sections (Prop. 12).

Because the bounds are also polar to each other (from Prop. 10), if we have a performance guarantee for $\Theta^{\circ}$, i.e.,

for some $\kappa\geqslant 1$, we immediately get

Proof We only consider the case where ${\mathcal{V}}$ is the unit disk, which gives the main idea of the proof of . We consider a maximizer $M$ of Eq. (40) and a vector $u=D^{1/2}\mathop{\rm sign}(v)$ where $v$ is a random sample from a normal distribution with mean and covariance matrix $D^{-1/2}MD^{-1/2}$, with $D=\mathop{\rm Diag}(\mathop{\rm diag}(M))$. We have, using standard arguments from , $uu^{\top}\in{\mathcal{C}}_{\mathcal{U}}\cap\{{\rm rank}(M)=1\}={\mathcal{P}}_{\mathcal{U}}$ and

The general case is also based on sampling and normalizing Gaussian random variables (see more details in ).

Assume that for all $w\in{\rm span}({\mathcal{U}})$ and $z\in{\rm span}({\mathcal{V}})$, we have the representations $\max\{\gamma_{\mathcal{U}}^{\circ}(w),\gamma_{\mathcal{U}}^{\circ}(-w)\}^{2}=\gamma_{{\mathcal{C}}_{\mathcal{U}}}^{\circ}(ww^{\top})$ and $\max\{\gamma_{\mathcal{V}}^{\circ}(z),\gamma_{\mathcal{V}}^{\circ}(-z)\}^{2}=\gamma_{{\mathcal{C}}_{\mathcal{V}}}^{\circ}(zz^{\top})$. Then Eq. (41) and Eq. (42) are valid for $\kappa=\min\{n,d\}$. If moreover, ${\mathcal{V}}$ is the unit disk, they are also valid for $\kappa=\sqrt{\min\{n,d\}}$.

Proof Our assumptions imply by duality that for all $w\in{\rm span}({\mathcal{U}})$, $\gamma_{{\mathcal{U}}\cup(-{\mathcal{U}})}(w)^{2}=\inf_{N\in{\mathcal{C}}_{\mathcal{U}}}w^{\top}N^{-1}w=G_{\mathcal{U}}(ww^{\top})$, with $G_{\mathcal{U}}(M)=\inf_{N\in{\mathcal{C}}_{\mathcal{U}}}MN^{-1}$, which is defined as, with $M=WW^{\top}\in{\mathcal{C}}_{\mathcal{U}}$, the minimum of $\mathop{\rm tr}B$ s.t. \Big{(}\!\begin{array}[]{cc}A&W\\ W^{\top}&N\end{array}\!\Big{)}\succcurlyeq 0. By taking $N=WW^{\top}$ and $A=I$, we obtain that $G_{\mathcal{U}}(M)\leqslant{\rm rank}(M)$.

Thus, there exists $u,v$ such that $\gamma_{\mathcal{U}}(u)>0$, $\gamma_{\mathcal{V}}(v)>0$ and u^{\top}Yv-\frac{\rho}{\min\{n,d\}}\big{(}\frac{1}{2}\gamma_{\mathcal{U}}(u)^{2}+\frac{1}{2}\gamma_{\mathcal{V}}(v)^{2}\big{)}. This implies that Eq. (41) is valid for $\kappa=\min\{n,d\}$.

Note that we may also obtain the same results without randomization. When $\gamma_{\mathcal{V}}=\|\cdot\|_{2}$, we have:

Thus by taking $v$ the largest eigenvector of $Y^{\top}MY$ and $u\in\mathop{\rm argmax}_{u\in{\mathcal{U}}}u^{\top}Yv$, we have

with an approximation ratio of $\min\{n,p\}$ between the operator norm and the nuclear norm.

The two previous propositions lead to candidates for vectors $u$ and $v$, through sampling from normal distributions with mean zero and covariance matrix equal to $M$ (or a normalized version thereof). Note that other possibilies are available, like taking a largest eigenvector (such as done in the proof of Prop. 12) and running the power method from Section 4.5 starting from any of the above candidates.

2 Random sampling

The following proposition provides an approximation ratio with high probability.

Note that the result above is rather weak as on top of the scaling in $\sqrt{n}$, the equivalence constants between $\Theta$ and $\|\cdot\|_{F}$ and between $\gamma_{\mathcal{U}}$ and $\|\cdot\|_{2}$ may be large as well (though finite because ${\mathcal{U}}^{\circ}$ is compact).

Upper bounds on approximation performance.

Obtaining decompositions from convex relaxations

The first possibility is to follow Section 4.4 and use solutions $Q,S$ of

and decompose $Q=AA^{\top}$ with $A$ full rank (i.e., $A^{\top}A$ invertible), and $S=BB^{\top}$ with $B$ full rank. Then $U=(AA^{\top})^{-1}AG\mathop{\rm Diag}(s)^{1/2}$ and $V=(BB^{\top})^{-1}BH\mathop{\rm Diag}(s)^{1/2}$, are such that $X=UV^{\top}$ and $\frac{1}{2}\gamma_{{\mathcal{C}}_{\mathcal{U}}}(UU^{\top})+\frac{1}{2}\gamma_{{\mathcal{C}}_{\mathcal{V}}}(VV^{\top})$ is minimal (and equal to the value of the relaxation). However, for any orthogonal matrix $R$, $(UR,VR)$ is also such a pair, and it is not possible to obtain a decomposition of $X$ such that the gauge functions $\gamma_{\mathcal{U}}$ and $\gamma_{\mathcal{V}}$ of all columns of $U$ and $V$ are small.

2 Conditional gradient algorithms

We may also find decompositions by approximately solving the following convex optimization problem (which is a generalized basis pursuit problem):

for $\lambda$ small enough, using convex optimization techniques that only access $\Theta$ through computing $\Theta^{\circ}(Y)=\sup_{(u,v)\in{\mathcal{U}}\times{\mathcal{V}}}u^{\top}Yv$, and the associated minimizers. This is exactly what generalized conditional gradient algorithms can do . However, we need an algorithm which is robust to obtaining only approximate maximizers $(u,v)$, with potentially multiplicative approximation guarantees for the computation of $\Theta^{\circ}$.

We consider the following algorithm started from $Z_{0}=0$, which iterates the following recursion, for $\rho_{t}=2/(t+1)$, $t\geqslant 1$:

In Appendix A, we show that if we can find only approximate maximizers $(u,v)$ with approximation ratio $\kappa\geqslant 1$, then, if $X_{\lambda}$ is the unique solution of Eq. (44), then we have

and $Z_{t}$ is a positive linear combination of matrices $u_{s-1}v_{s-1}^{\top}$, $s\leqslant t$, with a sum of coefficients which is less than $\Theta(Z_{t})$. The previous inequality implies that $\Theta(Z_{t})\leqslant\kappa\Theta(X)+O(1/(\lambda t)).$ Thus, when $\lambda$ is small enough and $t$ is large enough, we obtain an approximation of $\Theta(X)$ with approximation ratio which converges to a value less than $\kappa$.

Finding such decomposition by greedily and iteratively adding factors has been studied thoroughly in signal processing and statistics . In particular, if we assume that $\Theta$ is a norm, with our set of assumptions, the matching pursuit algorithm of may obtain an $\varepsilon$-approximation of $X$ with $O(\log\frac{1}{\varepsilon})$ rank-one factors, however, while the norm of the coefficients is bounded, it is not related to the decomposition gauge function $\Theta(X)$. In Appendix A.1, we show how optimizing over the scalar $\rho$ in the algorithm above leads to a similar result, while the sum of coefficients converge to a value which is less than $\kappa\Theta(X)$.

3 Simplicial methods

Conditional gradient algorithms to solve Eq. (44) may be extended by simply replacing step (b), which is the minimization over a half-line, by the minimization with respect to the cone generated by the already obtained rank-one matrices, i.e.,

This approach is sometimes referred to as fully corrective and is an instance of a simplicial method (see, e.g., ); typically, it requires much fewer iterations while the cost of each iteration is higher. Note that when the algorithm stops (i.e., there is no further progress in reducing the cost function), then we have solved Eq. (44) up to $\lambda(\kappa-1)\Omega(x_{\ast})$. An inbetween alternative is to optimize only over $\alpha$ and $\rho$. Note that the bound derived above also applies to these two extensions, which typically converge much quicker (see examples in Section 7).

Note that when $(u,v)$ is obtained from randomized rounding, the algorithm is related to what is proposed by , which uses non-adaptive random sampling for $u_{s-1}v_{s-1}^{\top}$. Moreover, the algorithm may be accelerated by only storing only vectors $u_{1},\dots,u_{t-1}$, and replacing the subproblem in Eq. (45) by

Simulations

We first compared several approaches to estimating $\Theta^{\circ}(Y)=\max_{u\in{\mathcal{U}},\ v\in{\mathcal{V}}}u^{\top}Yv=\max_{u\in\{0,1\}^{n}}\|Y^{\top}u\|_{2}$, for $Y$ a random matrix with independent and identically distributed components from a normal distribution with mean zero and variance one. We consider several strategies: (a) random sampling of $u$, then running the power method to convergence, (b) sampling from the solution of the relaxed semi-definite program, with and without running the power method, and (c) taking the non-randomized approach described in the proof of Prop. 12. In Figure 1, we can see that (a) the performance of the semidefinite-relaxation, even without the power method, is typically much better than the guarantee, (b) that sampling from the relaxed solution and then running the power method outperforms random initiatializations, and (c) the non-randomized rounding based on eigenvectors has a less stable behavior and sometimes performs better.

In Figure 1, we report averaged value of the randomized rounding procedures. If we take the best values over more than a thousand samples, the obtained values of $u^{\top}Yv$ of all three schemes happens to be very close, with a slight advantage to the initializations of the power methods from the convex relaxation.

Finding decomposition.

We aim to solve the problem in Eq. (44) for $\lambda=10^{-4}$, and consider four approaches: (a) using the semidefinite relaxation to obtain a lower-bound (with no explicit decomposition), (b) using the conditional gradient algorithm, (c) using a simplicial method (the regular version based on Eq. (45) or the one adapted to storing only values of $u$ in Eq. (46)), (d) using alternating optimization and (e) random sampling (Section 5.2). In Figure 3, we compare these algorithms in two situations, one where the relaxation is tight (left: $n=32$, $d=1$) and one where it is not tight (right: $n=32$, $d=16$). The simplicial algorithms are the faster to converge, with a clear advantage to the one that stores only the vectors $u$. The random selection procedure starts slow but eventually catches up, but never reaches the objective function of adaptive methods.

Finally, in Figure 2, we compare the result of using the simplicial method to obtain a decomposition to a simple alternating optimization method. We see that the convex relaxation outperforms significantly the non-convex approach.

Conclusion

We have presented a general framework for structured matrix decompositions based on gauge functions and semi-definite programming. Emphasis was put on situations where the rank-one factors belong to potentially non-convex and non-centrally symmetric sets. A series of algorithms and relaxations have been presented for a variety of structures.

Our limited experiments have focused on $\{0,1\}$-valued factors. It would be worth studying more precisely recovery guarantees and potentially more scalable algorithms for such cases, in particular for additional constraints such as cardinality or submodularity. Similarly, our general framework leads to altenatives to non-convex algorithms for non-negative matrix factorization, for which it may be possible to show tightness under certain assumptions similar to .

Moreover, it would be interesting to see if techniques designed to adaptively reduce the rank for nuclear-norm penalized problems may be extended to our setting as well in order to enforce a smaller number of factors. Finally, the non-convex approaches (power method and alternating minimization) tend to work well in practice, and finding sufficient conditions or random instances (see, e.g., for such work for sparse principal component analysis) where they provably behave well would provide valuable additional insights into the problem of structured matrix factorizations.

Acknowledgments

This work was supported by the European Research Council (SIERRA project 239993). The author thanks Guillaume Obozinski and Alexandre d’Aspremont for fruitful discussions related to this work.

Appendix A Generalized conditional gradient with approximate oracle

This implies by recursion (see for details) that

If moreover, $f$ is $\mu$-strongly convex and the line search is performed exactly (and without a bound on $\alpha$), then we may show a different bound. Up to the oracle with multiplicative approximation guarantees, the algorithm is then the same as what is proposed by , but our results use a slightly different set of assumptions.

First, we show that $f(x_{t})$ remains bounded. Indeed, we have $f(x_{t})\leqslant f((1-\rho_{t})x_{t-1})\leqslant(1-\rho_{t})f(x_{t-1})+\rho_{t}f(0)$, which leads to $f(x_{t})\leqslant f(0)$ for all $t\geqslant 1$. This implies that $f(0)\geqslant f(0)+x_{t}^{\top}f^{\prime}(0)+\frac{\mu}{2}\|x_{t}\|_{2}^{2}\geqslant-\|x_{t}\|_{2}\|f^{\prime}(0)\|_{2}+\frac{\mu}{2}\|x_{t}\|_{2}^{2}$, leading to $\|x_{t}\|_{2}\leqslant 2\|f^{\prime}(0)\|_{2}/\mu$.

We may then derive a different recursion, leading to

When applied to $f(x)=\frac{1}{2}\|x-y\|_{2}^{2}$, we obtain a decaying factor of

In this case, our algorithm is strongly related to the relaxed greedy algorithm of and our analysis provides an explicit link between these algorithms and basis pursuit .

We now assume that $f$ is $\mu$-strongly convex and that $f$ has a global minimum attained at $x_{\ast}$ such that $\gamma_{\mathcal{C}}(x_{\ast})\leqslant\omega$, and that we have a $\kappa$-approximate oracle for maximizing linear functions on ${\mathcal{C}}$. We consider the algorithm:

Moreover, we have $\displaystyle f(x_{t-1})-f(x_{\ast})\geqslant\frac{\mu}{2}\|x_{t-1}-x_{\ast}\|_{2}^{2}\geqslant\frac{\mu}{2}\frac{\gamma_{\mathcal{C}}(x_{t-1}-x_{\ast})^{2}}{\max_{\|u\|_{2}=1}\gamma_{\mathcal{C}}(u)^{2}}$, leading to, with $\Delta_{t}=f(x_{t-1})-f(x_{\ast})$,

If \displaystyle\bigg{[}\Delta_{t-1}+\frac{\varepsilon\sqrt{\mu}\Delta_{t-1}^{1/2}}{\sqrt{2}\max_{\|u\|_{2}=1}\gamma_{\mathcal{C}}(u)}\bigg{]}\frac{1}{L\kappa^{2}\omega^{2}\max_{u\in{\mathcal{C}}}\|u\|^{2}}<1, we have a minimizer $\rho\in[0,1)$, and

Thus, with \displaystyle\tau=\min\bigg{\{}\frac{1}{2},\frac{\varepsilon^{2}}{4\kappa^{2}\omega^{2}}\frac{\mu}{L}\frac{1}{\max_{\|u\|_{2}=1}\gamma_{\mathcal{C}}(u)\times\max_{u\in{\mathcal{C}}}\|u\|^{2}}\bigg{\}}, we have $\Delta_{t}\leqslant(1-\tau)\Delta_{t-1}$, and hence a linear convergence rate.

Appendix B Maximum of beta random variables

Using standard results from probability (see, e.g., ), we get, with $t=n/4$:

References