{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Minimizing Condition Number by Scaling\n",
    "This notebook provides and example of how to minimize the condition number of a matrix by scaling, taken from https://web.stanford.edu/~boyd/lmibook/. The problem is formulated in CVXPY as a generalized eigenvalue problem (GEVP) and utilizes the [DQCP](https://www.cvxpy.org/tutorial/dqcp/index.html) capabilities in CVXPY. The problem for a matrix $M \\in \\mathcal{R}^{m \\times n}$, with $m \\ge n$ is\n",
    "\n",
    "$$\n",
    "\\begin{array}{llr}\n",
    "\\text{minimize}   & \\gamma^2 & \\\\\n",
    "\\text{subject to} & P \\in \\mathcal{R}^{m \\times m} \\text{ and diagonal}, & P > 0, \\\\\n",
    "                  & Q \\in \\mathcal{R}^{n \\times n} \\text{ and diagonal}, & Q > 0, \\\\\n",
    "                  & Q \\le M^T P M \\le \\gamma^2 Q &\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example\n",
    "In the following code, we solve the above GEVP with CVXPY."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Is the problem DQCP?  True\n",
      "Initial Condition Number:  4.1538811703979786\n",
      "Optimized Condition Number:  1.7548711807791855\n"
     ]
    }
   ],
   "source": [
    "# import packages\n",
    "import cvxpy as cp\n",
    "import numpy as np\n",
    "\n",
    "# create helper functions\n",
    "def cond(A):\n",
    "    return np.linalg.cond(A)\n",
    "\n",
    "def evalCond(M,Q,P):\n",
    "    L     = np.diag(np.diag(P.value)**(1/2))\n",
    "    R     = np.diag(np.diag(Q.value)**(-1/2))\n",
    "    Mnew  = L @ M @ R\n",
    "    return np.linalg.cond(Mnew)\n",
    "\n",
    "# create a random matrix\n",
    "m = 3\n",
    "n = 2\n",
    "np.random.seed(2)\n",
    "M = np.random.rand(m,n)\n",
    "\n",
    "# specify the variables\n",
    "p = cp.Variable(m,pos=True)\n",
    "P = cp.diag(p)\n",
    "q = cp.Variable(n,pos=True)\n",
    "Q = cp.diag(q)\n",
    "\n",
    "# define the variables for GEVP\n",
    "A = M.T @ P @ M\n",
    "B = Q\n",
    "C = A - Q\n",
    "\n",
    "# create the constraints and objective\n",
    "ep = 1e-3\n",
    "constr = [C >= ep*np.eye(C.shape[0]),\n",
    "          P >= ep*np.eye(P.shape[0]),\n",
    "          Q >= ep*np.eye(Q.shape[0])]\n",
    "\n",
    "# note: the variable lambda = gamma^2 from the problem statement\n",
    "objFun = cp.Minimize(cp.gen_lambda_max(A,B))\n",
    "\n",
    "# create the problem\n",
    "problem = cp.Problem(objFun,constr)\n",
    "\n",
    "# check if DQCP\n",
    "print(\"Is the problem DQCP? \",problem.is_dqcp())\n",
    "\n",
    "# solve\n",
    "problem.solve(qcp=True,solver=\"SCS\")\n",
    "\n",
    "# print results\n",
    "if problem.status not in [\"infeasible\", \"unbounded\"]:\n",
    "    print(\"Initial Condition Number: \",cond(M))\n",
    "    print(\"Optimized Condition Number: \",evalCond(M,Q,P))\n",
    "else:\n",
    "    print(problem.status)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "cvxpy-dev",
   "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
