{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rank-one nonnegative matrix factorization\n",
    "\n",
    "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\n",
    "as the outer product of two positive vectors.\n",
    "\n",
    "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\n",
    "$xy^T$, that is,\n",
    "\n",
    "$$\n",
    "\\frac{1}{mn} \\sum_{i=1}^{m} \\sum_{j=1}^{n} R(A_{ij}, x_iy_j),\n",
    "$$\n",
    "\n",
    "where $R$ is the relative deviation of two positive numbers, defined as\n",
    "\n",
    "$$\n",
    "R(a, b) = \\max\\{a/b, b/a\\} - 1.\n",
    "$$\n",
    "\n",
    "The corresponding optimization problem is\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize} & \\frac{1}{mn} \\sum_{i=1}^{m} \\sum_{j=1}^{n} R(X_{ij}, x_iy_j)\n",
    "\\\\\n",
    "\\mbox{subject to} & x_1x_2 \\cdots x_m = 1 \\\\\n",
    "& X_{ij} = A_{ij}, \\quad \\text{for } (i, j) \\in \\Omega,\n",
    "\\end{array}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "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.\n",
    "\n",
    "The below code constructs and solves this optimization problem, with specific problem data\n",
    "\n",
    "$$\n",
    "A = \\begin{bmatrix}\n",
    "1.0 & ? &  1.9 \\\\\n",
    "? & 0.8 &  ? \\\\\n",
    "3.2 & 5.9&  ?\n",
    "\\end{bmatrix},\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal value:\n",
      " 1.7763568394002505e-14 \n",
      "\n",
      "Outer product approximation\n",
      " [[1.         1.84375    1.9       ]\n",
      " [0.43389831 0.8        0.82440678]\n",
      " [3.2        5.89999999 6.07999999]] \n",
      "\n",
      "x:  [0.89637009 0.38893346 2.86838428]\n",
      "y:  [1.11561063 2.0569071  2.1196602 ]\n"
     ]
    }
   ],
   "source": [
    "import cvxpy as cp\n",
    "\n",
    "m = 3\n",
    "n = 3\n",
    "X = cp.Variable((m, n), pos=True)\n",
    "x = cp.Variable((m,), pos=True)\n",
    "y = cp.Variable((n,), pos=True)\n",
    "\n",
    "outer_product = cp.vstack([x[i] * y for i in range(m)])\n",
    "relative_deviations = cp.maximum(\n",
    "  cp.multiply(X, outer_product ** -1),\n",
    "  cp.multiply(X ** -1, outer_product))\n",
    "objective = cp.sum(relative_deviations)\n",
    "constraints = [\n",
    "  X[0, 0] == 1.0,\n",
    "  X[0, 2] == 1.9,\n",
    "  X[1, 1] == 0.8,\n",
    "  X[2, 0] == 3.2,\n",
    "  X[2, 1] == 5.9,\n",
    "  x[0] * x[1] * x[2] == 1.0,\n",
    "]\n",
    "problem = cp.Problem(cp.Minimize(objective), constraints)\n",
    "problem.solve(gp=True)\n",
    "\n",
    "print(\"Optimal value:\\n\", 1.0/(m * n) * (problem.value - m * n), \"\\n\")\n",
    "print(\"Outer product approximation\\n\", outer_product.value, \"\\n\")\n",
    "print(\"x: \", x.value)\n",
    "print(\"y: \", y.value)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
