Rank-one nonnegative matrix factorization

The DGP atom library has several functions of positive matrices, including the trace, (matrix) product, sum, Perron-Frobenius eigenvalue, and \((I - X)^{-1}\) (eye-minus-inverse). In this notebook, we use some of these atoms to approximate a partially known elementwise positive matrix as the outer product of two positive vectors.

We would like to approximate \(A\) as the outer product of two positive vectors \(x\) and \(y\), with \(x\) normalized so that the product of its entries equals \(1\). Our criterion is the average relative deviation between the entries of \(A\) and \(xy^T\), that is,

\[\frac{1}{mn} \sum_{i=1}^{m} \sum_{j=1}^{n} R(A_{ij}, x_iy_j),\]

where \(R\) is the relative deviation of two positive numbers, defined as

\[R(a, b) = \max\{a/b, b/a\} - 1.\]

The corresponding optimization problem is

\[\begin{split}\begin{equation} \begin{array}{ll} \mbox{minimize} & \frac{1}{mn} \sum_{i=1}^{m} \sum_{j=1}^{n} R(X_{ij}, x_iy_j) \\ \mbox{subject to} & x_1x_2 \cdots x_m = 1 \\ & X_{ij} = A_{ij}, \quad \text{for } (i, j) \in \Omega, \end{array} \end{equation}\end{split}\]

with variables \(X \in \mathbf{R}^{m \times n}_{++}\), \(x \in \mathbf{R}^{m}_{++}\), and \(y \in \mathbf{R}^{n}_{++}\). We can cast this problem as an equivalent generalized geometric program by discarding the \(-1\) from the relative deviations.

The below code constructs and solves this optimization problem, with specific problem data

\[\begin{split}A = \begin{bmatrix} 1.0 & ? & 1.9 \\ ? & 0.8 & ? \\ 3.2 & 5.9& ? \end{bmatrix},\end{split}\]

import cvxpy as cp

m = 3
n = 3
X = cp.Variable((m, n), pos=True)
x = cp.Variable((m,), pos=True)
y = cp.Variable((n,), pos=True)

outer_product = cp.vstack([x[i] * y for i in range(m)])
relative_deviations = cp.maximum(
  cp.multiply(X, outer_product ** -1),
  cp.multiply(X ** -1, outer_product))
objective = cp.sum(relative_deviations)
constraints = [
  X[0, 0] == 1.0,
  X[0, 2] == 1.9,
  X[1, 1] == 0.8,
  X[2, 0] == 3.2,
  X[2, 1] == 5.9,
  x[0] * x[1] * x[2] == 1.0,
]
problem = cp.Problem(cp.Minimize(objective), constraints)
problem.solve(gp=True)

print("Optimal value:\n", 1.0/(m * n) * (problem.value - m * n), "\n")
print("Outer product approximation\n", outer_product.value, "\n")
print("x: ", x.value)
print("y: ", y.value)

Optimal value:
 1.7763568394002505e-14

Outer product approximation
 [[1.         1.84375    1.9       ]
 [0.43389831 0.8        0.82440678]
 [3.2        5.89999999 6.07999999]]

x:  [0.89637009 0.38893346 2.86838428]
y:  [1.11561063 2.0569071  2.1196602 ]