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 :math:(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 :math:A as the outer product of two positive vectors :math:x and :math:y, with :math:x normalized so that the product of its entries equals :math:1. Our criterion is the average relative deviation between the entries of :math:A and :math:xy^T, that is, .. math:: \frac{1}{mn} \sum_{i=1}^{m} \sum_{j=1}^{n} R(A_{ij}, x_iy_j), where :math:R is the relative deviation of two positive numbers, defined as .. math:: R(a, b) = \max\{a/b, b/a\} - 1. The corresponding optimization problem is .. math:: $$\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}$$ with variables :math:X \in \mathbf{R}^{m \times n}_{++}, :math:x \in \mathbf{R}^{m}_{++}, and :math:y \in \mathbf{R}^{n}_{++}. We can cast this problem as an equivalent generalized geometric program by discarding the :math:-1 from the relative deviations. The below code constructs and solves this optimization problem, with specific problem data .. math:: A = \begin{bmatrix} 1.0 & ? & 1.9 \\ ? & 0.8 & ? \\ 3.2 & 5.9& ? \end{bmatrix}, .. code:: python 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) .. parsed-literal:: 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 ]