Matrix factorization is a fundamental tool in data mining and machine learning, with applications in clustering, topic modeling, recommender systems, and interpretable representation learning. Given a data matrix $\mathbf{X}\in\mathbb{R}^{d\times n}$, a typical factorization model seeks a low-rank approximation

where $\fm\in\mathbb{R}^{d\times k}$ contains basis components and $\gm\in\mathbb{R}^{k\times n}$ contains low-dimensional representations of the samples. Classical nonnegative matrix factorization (NMF) improves interpretability by imposing nonnegativity on both factors, while many variants further introduce sparsity, orthogonality, or graph regularization.

Existing factorization methods still face two limitations. First, the learned basis components can be highly redundant. In NMF and related models, different columns of $\fm$ may become strongly correlated, making the learned representation less interpretable and less discriminative for downstream tasks such as clustering. Second, many classical algorithms, especially multiplicative-update methods, are tightly coupled with nonnegativity assumptions and are therefore less flexible for real-valued data or more general constrained factorization settings. These observations motivate a matrix factorization framework that learns diverse basis components, produces interpretable sample representations, and remains applicable beyond standard NMF.

Simplex-Constrained Incoherent Matrix Factorization

In this blog, we will discuss Simplex-Constrained Incoherent Matrix Factorization. The key idea is to represent each sample as a convex combination of basis vectors by constraining every column of $\gm$ to lie on the probability simplex:

Thus, $ \vx_j\approx \fm\vg_j=\sum_{r=1}^k g_{rj}\vf_r, $ where $g_{rj}$ measures the contribution of the $r$-th basis component to the $j$-th sample. This simplex constraint gives $\gm$ a natural probabilistic interpretation as a soft cluster-assignment matrix and makes the representation useful for clustering and interpretable data analysis. In contrast, K-means assigns each sample to a hard one-hot indicator vector. Thus, our formulation can be viewed as a simplex-based soft relaxation of the hard assignment used in K-means.

In addition, to encourage diverse basis components, we introduce an incoherence regularizer based on the positive off-diagonal entries of the Gram matrix:

Therefore, the regularizer discourages positive correlation among different basis vectors while not penalizing negative inner products. Compared with a linear inner-product penalty, the squared positive-part formulation places stronger emphasis on highly positively correlated basis pairs and provides a simple smooth penalty away from the hinge point.

We further impose a data-adaptive column norm constraint on $\fm$: $ |\vf_i|_2\le R, \ i=1,\dots,k, $ where $R$ can be chosen according to the scale of the data, such as $R=\max_j|\vx_j|_2$. This constraint controls the scale of the basis matrix without requiring data normalization. When nonnegative bases are desired, we additionally impose $\fm\ge 0$; for real-valued data, we simply remove this nonnegativity constraint. Therefore, the same model naturally covers both NMF-style and real-valued matrix factorization settings.

The proposed optimization problem is

subject to

where $\mathcal{C}_F$ denotes either the nonnegative scale-constrained set or the real-valued scale-constrained set. To optimize this model, we develop a hybrid projected-gradient and entropy mirror-descent algorithm as belows.

Algorithm: Hybrid Projected–Entropic Mirror Descent

Input: Data matrix $\xm$, rank $k$, parameters $\lambda, R, \epsilon$

Output: Factor matrices $\fm, \gm$

1. Initialize $\fm^0 \in \mathcal{C}_F$ and $\vg_j^0 \in \operatorname{relint}(\Delta_k)$ for $j = 1,\ldots,n$.

2. For $t = 0,1,2,\ldots$:

   1. Compute

      \[\sm^t = \left[ (\fm^t)^\top \fm^t - \operatorname{diag}\big((\fm^t)^\top \fm^t\big) \right]_+ .\]

   2. Compute

      \[\nabla_{\fm}\Phi(\fm^t,\gm^t) = (\fm^t\gm^t-\xm)(\gm^t)^\top + 4\lambda \fm^t\sm^t .\]

   3. Set $\alpha_t = 1/\overline{L}_F$ and update

      \[\fm^{t+1} = \Pi_{\mathcal{C}_F} \left( \fm^t-\alpha_t\nabla_{\fm}\Phi(\fm^t,\gm^t) \right).\]

   4. Set

      \[\eta_t = \frac{1}{\|\fm^{t+1}\|_2^2+\epsilon}.\]

   5. For $j = 1,\ldots,n$:

      1. Compute

         \[\vq_j^{t+1} = (\fm^{t+1})^\top (\fm^{t+1}\vg_j^t-\vx_j).\]

      2. For $r = 1,\ldots,k$, update

         \[g_{rj}^{t+1} = \frac{ g_{rj}^{t}\exp(-\eta_t q_{rj}^{t+1}) }{ \sum_{\ell=1}^{k} g_{\ell j}^{t}\exp(-\eta_t q_{\ell j}^{t+1}) }.\]

   6. If a stopping criterion is satisfied, break.

3. Return $\fm^{t+1}, \gm^{t+1}$.

The proposed algorithm alternates between a projected-gradient step for $\fm$ and an entropy mirror-descent step for $\gm$. The $\fm$-update is Euclidean because $\mathcal{C}_F$ is a simple column-wise norm constraint. The $\gm$-update is entropic because each column $\vg_j$ lies on the probability simplex.

For the entropy mirror step, we use the negative entropy function

defined on the simplex. Its associated Bregman divergence is the Kullback–Leibler divergence

For the $\fm$-update, we use the Euclidean mirror map

whose associated Bregman divergence is

At iteration $t$, given $(\fm^t,\gm^t)$, we first update $\fm$ by

where $\Pi_{\mathcal{C}_F}$ is the Euclidean projection onto $\mathcal{C}_F$ which is separable across columns.

For the real-valued scale-constrained case,

the projection is

For the nonnegative case,

we first threshold negative entries and then project onto the $\ell_2$ ball:

For each column $\vg_j$, we first compute

and update $\vg_j$ by entropy mirror descent:

This subproblem has the closed-form multiplicative-normalization update [1]:

Thus, simplex feasibility is preserved automatically.

Convergence Analysis

We now introduce Pinsker’s inequality [2], which will be used in the convergence analysis. For any $\vu,\vv\in\Delta_k$,

Lemma 1. Fix $\fm$ and let $\vg^+$ be generated from $\vg \in \Delta_k$ by

If $0 < \eta < 1/L_G(\fm)$, then

Consequently,

Proof.
By the optimality of the mirror-descent subproblem, comparing the objective value at $\vg^+$ and $\vg$ gives

Hence,

Using the smoothness inequality with $\vu=\vg^+$ and $\vv=\vg$, we obtain

By Pinsker’s inequality,

Therefore,

which proves the first part. Applying Pinsker’s inequality once more gives the second part.

Corollary 1: Fix $\fm$ and update all columns of $\gm$, if $0<\eta<1/|\fm|_2^2$, then

Consequently,

Lemma 2: Let

If $0<\alpha<2/L_F(\gm)$, then [3]

Theorem 1: Suppose that the stepsizes satisfy

Then the sequence generated by Algorithm 1 satisfies

More precisely,

where

Lemma 3: A uniform upper bound for $\fm$ update is

Theorem 2: Assume that the stepsizes are chosen such that there exist constants $c_F>0$ and $c_G>0$ satisfying

for all $t$. Then the sequence ${(\fm^t,\gm^t)}$ generated by Algorithm 1 is bounded, the objective values ${\Phi(\fm^t,\gm^t)}$ converge, and

Moreover, every accumulation point $(\fm^\star,\gm^\star)$ is a stationary point in the sense that

and

Experimental Results

References

[1] Nisheeth Vishnoi. “Algorithms for convex optimization.” Cambridge University Press, 2021.

[2] Clément L. Canonne. “A short note on an inequality between KL and TV”.

[3] Amir Beck. “First-order methods in optimization”.