{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PCA\n",
    "\n",
    "In this notebook, I will review a little bit about [PCA](https://arxiv.org/abs/1404.1100), implement [recursive PCA](http://www.sciencedirect.com/science/article/pii/S0959152400000226), how PCA can be viewed as an optimization problem, and implement a constrained version of this optimization process for PCA applied on time-series by including a roughness penalty.\n",
    "\n",
    "[Generalized PCA](https://arxiv.org/abs/1202.4002) will not be considered after careful consideration. Note: \"Generalized PCA\" is \"an algebro-geometric solution to the problem of segmenting an unknown number of subspaces of unknown and varying dimensions from sample data points\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib.patches import FancyArrowPatch\n",
    "from mpl_toolkits.mplot3d import proj3d\n",
    "\n",
    "from sklearn.decomposition import PCA\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Arrow3D(FancyArrowPatch):\n",
    "    def __init__(self, xs, ys, zs, *args, **kwargs):\n",
    "        FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)\n",
    "        self._verts3d = xs, ys, zs\n",
    "\n",
    "    def draw(self, renderer):\n",
    "        xs3d, ys3d, zs3d = self._verts3d\n",
    "        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)\n",
    "        self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))\n",
    "        FancyArrowPatch.draw(self, renderer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first start to apply PCA on the following toy example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Toy data\n",
    "X = np.array([[-1, -1],\n",
    "              [-2, -1],\n",
    "              [-3, -2],\n",
    "              [1, 1],\n",
    "              [2, 1],\n",
    "              [3, 2]], dtype=np.float64) # NxT\n",
    "\n",
    "# Plot\n",
    "plt.figure(figsize=(5,4))\n",
    "plt.scatter(X[:,0], X[:,1], color='b')\n",
    "plt.arrow(0, 0, 1, 0, length_includes_head = True, head_width = 0.15, color='k')\n",
    "plt.arrow(0, 0, 0, 1, length_includes_head = True, head_width = 0.15, color='')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### PCA\n",
    "n_components = 2 # number of principal axes that we want to keep\n",
    "\n",
    "## Let's compute PCA from scratch\n",
    "# 1. Center the data\n",
    "mean = X.mean(axis=0)\n",
    "X -= mean # note that the data was already centered\n",
    "N = X.shape[0]\n",
    "\n",
    "# 2. Compute the covariance matrix\n",
    "CovX = 1./(N-1) * X.T.dot(X) # TxT (same as np.cov(X, rowvar=False)))\n",
    "\n",
    "# 3. Compute the eigenvectors of this covariance matrix\n",
    "# np.linalg.eigh is more efficient than np.linalg.eig for symmetric matrix\n",
    "evals, evecs = np.linalg.eigh(CovX)\n",
    "\n",
    "# 4. Sort the eigenvalues (in decreasing order) and eigenvectors\n",
    "idx = np.argsort(evals)[::-1]\n",
    "evals = evals[idx]\n",
    "evecs = evecs[:,idx]\n",
    "\n",
    "# 5. Form the projection matrix\n",
    "P = evecs[:,:n_components]\n",
    "print(P)\n",
    "\n",
    "# 5. Project the data using the projection matrix\n",
    "# This is the same as rotating the matrix X using P\n",
    "Y = X.dot(P)\n",
    "    \n",
    "# 6. Compare it with standard PCA\n",
    "pca = PCA(n_components=n_components)\n",
    "pca = pca.fit(X)\n",
    "Xnew = pca.transform(X)\n",
    "\n",
    "print(pca.components_.T)\n",
    "print(np.allclose(Xnew, Y))\n",
    "\n",
    "# 7. Plot the data\n",
    "plt.figure(figsize=(10,4))\n",
    "plt.subplot(1,2,1)\n",
    "plt.title('PCA: eigenvalues')\n",
    "plt.bar(np.array([0.,0.1]), evals, width=0.1)\n",
    "plt.xlim(0.,1.)\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.title('PCA: data and eigenvectors')\n",
    "plt.scatter(X[:,0], X[:,1], color='b')\n",
    "plt.arrow(0, 0, np.sqrt(evals[0])*P[0,0], np.sqrt(evals[0])*P[1,0],\n",
    "          length_includes_head = True, head_width = 0.15, color='b')\n",
    "plt.arrow(0, 0, np.sqrt(evals[1])*P[0,1], np.sqrt(evals[1])*P[1,1],\n",
    "          length_includes_head = True, head_width = 0.15, color='b')\n",
    "plt.scatter(Y[:,0], Y[:,1], color='r')\n",
    "plt.arrow(0, 0, np.sqrt(evals[0]), 0, length_includes_head = True, head_width = 0.15, color='r')\n",
    "plt.arrow(0, 0, 0, np.sqrt(evals[1]), length_includes_head = True, head_width = 0.15, color='r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The covariance matrix can be recovered from these eigenvalues and eigenvectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the eigenvectors and eigenvalues, we can of course recover the covariance matrix\n",
    "# Note: if P is not square, we have to fill it up.\n",
    "np.allclose(CovX, P.dot(np.diag(evals)).dot(P.T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Few properties:\n",
    "\n",
    "* Applying PCA several times is the same as applying it one time. This is because PCA diagonalizes our matrix, thus applying PCA on a diagonal matrix will result in the same matrix.\n",
    "* Applying PCA on a part of the data and another PCA on the other part, then applying PCA on the concatenation of both do not result in the same matrix as applying PCA on the whole data.\n",
    "* Applying PCA on a rotated matrix does not give the same result as applying PCA on the initial matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here is the general method\n",
    "\n",
    "def pca(X, normalize=False, copy=True):\n",
    "    if copy:\n",
    "        X = np.copy(X)\n",
    "\n",
    "    # 1. Center the data\n",
    "    mean = X.mean(axis=0)\n",
    "    X -= mean\n",
    "    N = X.shape[0]\n",
    "    \n",
    "    if normalize:\n",
    "        X /= X.std(axis=0)\n",
    "\n",
    "    # 2. Compute the covariance matrix\n",
    "    CovX = 1./(N-1) * X.T.dot(X) # TxT (same as np.cov(X, rowvar=False)))\n",
    "\n",
    "    # 3. Compute the eigenvectors of this covariance matrix\n",
    "    # np.linalg.eigh is more efficient than np.linalg.eig for symmetric matrix\n",
    "    evals, evecs = np.linalg.eigh(CovX)\n",
    "\n",
    "    # 4. Sort the eigenvalues (in decreasing order) and eigenvectors\n",
    "    idx = np.argsort(evals)[::-1]\n",
    "    evals, evecs = evals[idx], evecs[:,idx]\n",
    "\n",
    "    return evals, evecs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Applying PCA on time series\n",
    "\n",
    "A fundamental question when applying PCA on time series is how to visualize this high dimensional data. Indeed, a sample $\\pmb{x}(t) \\in \\mathbb{R}^T$. One way is to plot this $\\pmb{x}(t)$ where the x-axis is the time, and y-axis is $x(t)$. Each time $t_i$ $(\\forall i \\in \\{0,...,T\\})$ represents a dimension. By plotting a vertical line at time $t=t_i$, we can see the variance in this particular dimension.\n",
    "\n",
    "For an infinite vector, or function, check about *Functional PCA*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recursive PCA\n",
    "\n",
    "Let's now apply **recursive PCA** on this toy example, with 3 new samples coming at different time steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's augment our matrix with 3 new samples\n",
    "Xs = np.array([[-1,1],\n",
    "               [-3,0],\n",
    "               [-4,-1]], dtype=np.float64)\n",
    "\n",
    "n_components = 2\n",
    "k = float(N)\n",
    "R = CovX\n",
    "X_aug = X\n",
    "mean = X.mean(axis=0).reshape(-1,1)\n",
    "print(evals)\n",
    "for x in Xs:\n",
    "    x = x.reshape(-1,1)\n",
    "    X_aug = np.vstack((X_aug, x.T))\n",
    "    X_aug1 = X_aug - X_aug.mean(axis=0)\n",
    "    pca = PCA(n_components=n_components)\n",
    "    pca = pca.fit(X_aug1)\n",
    "    #print(pca.get_covariance())\n",
    "\n",
    "    # Recursive PCA\n",
    "    new_mean = k/(k+1) * mean + 1./(k+1) * x\n",
    "    diff_mean = (new_mean - mean)\n",
    "    R = (k-1)/k * R + diff_mean.dot(diff_mean.T) + 1./k * (x-new_mean).dot((x-new_mean).T)\n",
    "    #print(R)\n",
    "    print(\"The cov of the whole augmented matrix is equal to the recursive cov: {0}\".format(\n",
    "            np.allclose(pca.get_covariance(), R)))\n",
    "    k+=1\n",
    "    mean = new_mean\n",
    "    \n",
    "    evals = np.linalg.eigh(R)[0]\n",
    "    idx = np.argsort(evals)[::-1]\n",
    "    evals = evals[idx]\n",
    "    print(evals)\n",
    "\n",
    "# Compute the new projection matrix\n",
    "evals, evecs = np.linalg.eigh(R)\n",
    "idx = np.argsort(evals)[::-1]\n",
    "evals, evecs = evals[idx], evecs[:,idx]\n",
    "P = evecs[:,:n_components]\n",
    "Y = X_aug.dot(P)\n",
    "    \n",
    "# Plot the data\n",
    "plt.figure(figsize=(10,4))\n",
    "plt.subplot(1,2,1)\n",
    "plt.title('PCA: eigenvalues')\n",
    "plt.bar(np.array([0.,0.1]), evals, width=0.1)\n",
    "plt.xlim(0.,1.)\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.title('PCA: data and eigenvectors')\n",
    "plt.scatter(X_aug[:,0], X_aug[:,1], color='b')\n",
    "plt.arrow(0, 0, np.sqrt(evals[0])*P[0,0], np.sqrt(evals[0])*P[1,0],\n",
    "          length_includes_head = True, head_width = 0.15, color='b')\n",
    "plt.arrow(0, 0, np.sqrt(evals[1])*P[0,1], np.sqrt(evals[1])*P[1,1],\n",
    "          length_includes_head = True, head_width = 0.15, color='b')\n",
    "plt.scatter(Y[:,0], Y[:,1], color='r')\n",
    "plt.arrow(0, 0, np.sqrt(evals[0]), 0, length_includes_head = True, head_width = 0.15, color='r')\n",
    "plt.arrow(0, 0, 0, np.sqrt(evals[1]), length_includes_head = True, head_width = 0.15, color='r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now add a sample that can not be modeled by a linear combination of the principal axes, i.e. which is orthogonal to the current covariance matrix. Then, as usual, let's apply recursive PCA on it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's add a sample that can not be modeled by a linear combination of the principal axes\n",
    "# i.e. which is orthogonal to the current covariance matrix.\n",
    "# Then, let's apply recursive PCA on it.\n",
    "\n",
    "# New 3D sample\n",
    "x = np.array([1,-1,1], dtype=np.float64).reshape(-1,1)\n",
    "\n",
    "# Reshaping previous values (pad a column/row of zeros)\n",
    "X_aug = np.pad(X_aug, ((0,0), (0,1)), mode='constant') # Nx(T+1)\n",
    "mean = np.pad(mean, ((0,1),(0,0)), mode='constant')\n",
    "R = np.pad(R, ((0,1),(0,1)), mode='constant')\n",
    "\n",
    "# Adding new sample and compute mean\n",
    "X_aug = np.vstack((X_aug, x.T))\n",
    "X_aug1 = X_aug - X_aug.mean(axis=0)\n",
    "\n",
    "# Use sklearn PCA\n",
    "n_components = 3\n",
    "pca = PCA(n_components=n_components)\n",
    "pca = pca.fit(X_aug1)\n",
    "#print(pca.get_covariance())\n",
    "\n",
    "# Recursive PCA\n",
    "new_mean = k/(k+1) * mean + 1./(k+1) * x\n",
    "diff_mean = (new_mean - mean)\n",
    "R = (k-1)/k * R + diff_mean.dot(diff_mean.T) + 1./k * (x-new_mean).dot((x-new_mean).T)\n",
    "#print(R)\n",
    "print(\"The cov of the whole augmented matrix is equal to the recursive cov: {0}\".format(\n",
    "            np.allclose(pca.get_covariance(), R)))\n",
    "k+=1\n",
    "mean = new_mean\n",
    "#print('-'*30)\n",
    "\n",
    "# Compute the new projection matrix\n",
    "evals, evecs = np.linalg.eigh(R)\n",
    "idx = np.argsort(evals)[::-1]\n",
    "evals, evecs = evals[idx], evecs[:,idx]\n",
    "P = evecs[:,:n_components]\n",
    "Y = X_aug.dot(P)\n",
    "    \n",
    "# Plot the data\n",
    "fig = plt.figure(figsize=(10,4))\n",
    "plt.subplot(1,2,1)\n",
    "plt.title('PCA: eigenvalues')\n",
    "plt.bar(np.array([0.,0.1,0.2]), evals, width=0.1)\n",
    "plt.xlim(0.,1.)\n",
    "\n",
    "ax = fig.add_subplot(122, projection='3d')\n",
    "ax.set_title('PCA: data and eigenvectors')\n",
    "ax.scatter(X_aug[:,0], X_aug[:,1], X_aug[:,2])\n",
    "# From https://stackoverflow.com/questions/22867620/putting-arrowheads-on-vectors-in-matplotlibs-3d-plot\n",
    "for v in evecs:\n",
    "    a = Arrow3D([0., v[0]], [0., v[1]], [0., v[2]],\n",
    "                mutation_scale=20, lw=1, arrowstyle=\"-|>\", color=\"b\")\n",
    "    ax.add_artist(a)\n",
    "ax.scatter(Y[:,0], Y[:,1], Y[:,2], color='r')\n",
    "a = Arrow3D([0., evals[0]], [0., evals[1]], [0., evals[2]],\n",
    "            mutation_scale=20, lw=1, arrowstyle=\"-|>\", color=\"r\")\n",
    "ax.add_artist(a)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## PCA as an Optimization Problem\n",
    "\n",
    "PCA can be viewed as an optimization problem in 2 different ways. Theses 2 approaches are equivalent.\n",
    "1. Maximize the variance of the projected data.\n",
    "2. Minimize the reconstruction error in a least-square sense.\n",
    "\n",
    "Mathematically, here is the optimization problem that we are trying to solve:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\max_{\\pmb{v_i}} \\: \\pmb{v_i}^T \\pmb{X}^T \\pmb{X v_i} \\quad \\mbox{subj. to} \\quad \\begin{array}{l} \\pmb{v_i}^T \\pmb{v_i} = 1 \\\\ \\pmb{v_i}^T \\pmb{v_j} = 0\n",
    "\\end{array},\n",
    "\\end{equation}\n",
    "$\\forall i \\in \\{1,...,D\\}, \\forall 1\\leq j < i$.\n",
    "\n",
    "Nice references:\n",
    "* [What is the objective fct of PCA? (StackExchange)](https://stats.stackexchange.com/questions/10251/what-is-the-objective-function-of-pca)\n",
    "* [PCA objective function: what is the connection between maximizing variance and minimizing error? (StackExchange)](https://stats.stackexchange.com/questions/32174/pca-objective-function-what-is-the-connection-between-maximizing-variance-and-m)\n",
    "* [Everything you did and didn't know about PCA (blog)](http://alexhwilliams.info/itsneuronalblog/2016/03/27/pca/)\n",
    "* [\"PCA and Optimization - A Tutorial\", 2015 (paper)](http://scholarscompass.vcu.edu/cgi/viewcontent.cgi?article=1006&context=ssor_pubs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data\n",
    "# Toy data\n",
    "X = np.array([[-1, -1],\n",
    "              [-2, -1],\n",
    "              [-3, -2],\n",
    "              [1, 1],\n",
    "              [2, 1],\n",
    "              [3, 2]], dtype=np.float64) # NxT\n",
    "N = X.shape[0]\n",
    "\n",
    "# PCA\n",
    "evals, evecs = pca(X)\n",
    "print(evals)\n",
    "print(evecs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using Scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize\n",
    "\n",
    "# PCA as an optimization\n",
    "\n",
    "# cache the 'covariance' matrix\n",
    "C = X.T.dot(X)/(N-1)\n",
    "\n",
    "# define objective function to MINIMIZE\n",
    "f = lambda u: -(u.T.dot(C)).dot(u)\n",
    "\n",
    "# define initial guess\n",
    "x0 = np.zeros((2,1))\n",
    "\n",
    "# define optimization method\n",
    "# By default, it will be 'BFGS', 'L-BFGS-B', or 'SLSQP' depending on the constraints and bounds\n",
    "method = None\n",
    "\n",
    "# define constraints\n",
    "constraints = [{'type': 'eq', 'fun': lambda u: u.T.dot(u) - 1}]\n",
    "\n",
    "# Minimize --> it returns an instance of OptimizeResult\n",
    "u1 = minimize(f, x0, method=method, constraints=constraints)\n",
    "#print(u1)\n",
    "u1 = u1.x.reshape(-1,1) # get solution\n",
    "\n",
    "# Add constraint\n",
    "constraints.append({'type': 'eq', 'fun': lambda u: u1.T.dot(u)})\n",
    "u2 = minimize(f, x0, method=method, constraints=constraints)\n",
    "#print(u2)\n",
    "u2 = u2.x.reshape(-1,1)\n",
    "\n",
    "# stack the optimized vector found\n",
    "P = np.hstack((u1,u2))\n",
    "print(P)\n",
    "\n",
    "# Plot\n",
    "plt.title('PCA: data and eigenvectors')\n",
    "plt.scatter(X[:,0], X[:,1], color='b')\n",
    "plt.arrow(0, 0, P[0,0], P[1,0], length_includes_head = True, head_width = 0.15, color='b')\n",
    "plt.arrow(0, 0, P[0,1], P[1,1], length_includes_head = True, head_width = 0.15, color='g')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define PCA optimization method\n",
    "def pca_scipy(X, rough_param=0.0, normalize=False, copy=True):\n",
    "    \"\"\"\n",
    "    Compute PCA on the given data using an optimization process.\n",
    "    \"\"\"\n",
    "    if copy:\n",
    "        X = np.copy(X)\n",
    "\n",
    "    # center the data\n",
    "    mean = X.mean(axis=0)\n",
    "    X -= mean\n",
    "    N = X.shape[0]\n",
    "    T = X.shape[1]\n",
    "\n",
    "    # normalize\n",
    "    if normalize:\n",
    "        X /= X.std(axis=0)\n",
    "\n",
    "    \n",
    "    # compute 'covariance' matrix and cache it\n",
    "    N = X.shape[0]\n",
    "    C = X.T.dot(X)/(N-1)\n",
    "\n",
    "    # define objective function to MINIMIZE\n",
    "    #f = lambda u: -(u.T.dot(C)).dot(u)\n",
    "    def f(u):\n",
    "        pen = 0\n",
    "        if rough_param != 0 and u.size > 2:\n",
    "            ddu = np.diff(np.diff(u))\n",
    "            rough_pen = (ddu**2).sum()\n",
    "            pen = rough_param*rough_pen\n",
    "            #if u.size > 4:\n",
    "            #    ddddu = np.diff(np.diff(ddu))\n",
    "            #    rough_pen = (ddddu**2).sum()\n",
    "            #    pen += rough_param*rough_pen\n",
    "        loss = -(u.T.dot(C)).dot(u)\n",
    "        return loss + pen\n",
    "\n",
    "    # define initial guess\n",
    "    x0 = np.ones((T,)) #np.zeros((T,))\n",
    "\n",
    "    # define optimization method\n",
    "    # By default, it will be 'BFGS', 'L-BFGS-B', or 'SLSQP' depending on the constraints and bounds\n",
    "    # If constraints, it will be 'SLSQP' (Sequential Least SQuares Programming)\n",
    "    # 'Nelder-Mead', 'Powell', 'CG', 'BFGS', 'Newton-CG', 'L-BFGS-B', 'TNC', 'COBYLA', 'SLSQP', 'dogleg', 'trust-ncg'\n",
    "    # 'Nelder-Mead', 'Powell', 'CG', 'Newton-CG', 'TNC', 'COBYLA', 'dogleg', 'trust-ncg' can not handle (eq) constraints\n",
    "    # 'BFGS', 'L-BFGS-B' do not work\n",
    "    method = 'SLSQP'\n",
    "\n",
    "    # define 1st constraints: norm of 1\n",
    "    constraints = [{'type': 'eq', 'fun': lambda u: u.T.dot(u) - 1}]\n",
    "\n",
    "    # optimize recursively\n",
    "    evals, evecs = [], []\n",
    "    messages = {}\n",
    "    for i in range(T):\n",
    "        if i != 0:\n",
    "            # add orthogonality constraint\n",
    "            constraints.append({'type': 'eq', 'fun': lambda u: u1.T.dot(u)})\n",
    "\n",
    "        # minimize --> it returns an instance of OptimizeResult\n",
    "        u1 = minimize(f, x0, method=method, constraints=constraints)\n",
    "        if not u1.success:\n",
    "            messages[i] = u1.message\n",
    "\n",
    "        # get 'eigenvalue'\n",
    "        evals.append(-u1.fun)\n",
    "\n",
    "        # get solution\n",
    "        u1 = u1.x\n",
    "        evecs.append(u1)\n",
    "\n",
    "    return np.array(evals), np.array(evecs).T, messages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "evals, evecs, messages = pca_scipy(X)\n",
    "P = evecs\n",
    "print(messages)\n",
    "print(evals)\n",
    "print(P)\n",
    "\n",
    "# Plot\n",
    "plt.figure(figsize=(10,4))\n",
    "plt.subplot(1,2,1)\n",
    "plt.title('PCA: eigenvalues')\n",
    "plt.bar(np.array([0.,0.1]), evals, width=0.1)\n",
    "plt.xlim(0.,1.)\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.title('PCA: data and eigenvectors')\n",
    "plt.scatter(X[:,0], X[:,1], color='b')\n",
    "plt.arrow(0, 0, P[0,0], P[1,0], length_includes_head = True, head_width = 0.15, color='b')\n",
    "plt.arrow(0, 0, P[0,1], P[1,1], length_includes_head = True, head_width = 0.15, color='g')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using CVXPY\n",
    "\n",
    "Note: You **cannot** use `cvxpy` for this problem, as we are trying to maximize a convex function, and `cvxpy` only accepts to maximize a concave fct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxpy as cvx\n",
    "\n",
    "# cache the 'covariance' matrix\n",
    "C = X.T.dot(X)/(N-1)\n",
    "\n",
    "# define vector to optimize\n",
    "u1 = cvx.Variable(X.shape[1])\n",
    "print(cvx.quad_form(u1, C).is_dcp())\n",
    "print(cvx.quad_form(u1, C).is_quadratic())\n",
    "\n",
    "# define objective fct to maximize\n",
    "#f = cvx.Maximize(u1.T*C*u1) \n",
    "f = cvx.Maximize(cvx.quad_form(u1, C)) # this does not work!\n",
    "#f = cvx.Minimize(cvx.quad_form(u1, C)) # this works (if no constraints) but that is not what we want to achieve!\n",
    "constraints = [u1.T*u1 == 1]\n",
    "prob = cvx.Problem(f, constraints)\n",
    "\n",
    "result = prob.solve()\n",
    "print(u1.value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using NLopt\n",
    "\n",
    "Nonlinear optimization algorithms that can handle nonlinear inequality and **equality** constraints are:\n",
    "* ISRES (Improved Stochastic Ranking Evolution Strategy) $\\rightarrow$ global derivative-free\n",
    "* COBYLA (Constrained Optimization BY Linear Approximations) $\\rightarrow$ local derivative-free\n",
    "* SLSQP (Sequential Least-SQuares Programming) $\\rightarrow$ local gradient-based\n",
    "* AUGLAG (AUGmented LAGrangian) $\\rightarrow$ global/local derivative-free/gradient based (determined based on the subsidiary algo)\n",
    "\n",
    "More information about the various algorithms can be found on this [link](https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Playground with NLopt\n",
    "import nlopt\n",
    "\n",
    "nlopt_results = {1: 'success', 2: 'stop_val reached', 3: 'ftol reached', 4: 'xtol reached',\n",
    "                 5: 'maxeval reached', 6: 'maxtime reached', -1: 'failure', -2: 'invalid args',\n",
    "                 -3: 'out of memory', -4: 'roundoff limited', -5: 'forced stop'}\n",
    "n_iter = 0\n",
    "N = X.shape[0]\n",
    "\n",
    "# cache the 'covariance' matrix\n",
    "C = X.T.dot(X)/(N-1)\n",
    "\n",
    "# define random seed\n",
    "nlopt.srand(125)\n",
    "\n",
    "# define which solver to use\n",
    "#optimizer = nlopt.GN_ISRES\n",
    "#optimizer = nlopt.LN_COBYLA\n",
    "optimizer = nlopt.LD_SLSQP\n",
    "#optimizer = nlopt.LD_AUGLAG # nlopt.AUGLAG, nlopt.AUGLAG_EQ, nlopt.LD_AUGLAG,\n",
    "                            # nlopt.LD_AUGLAG_EQ, nlopt.LN_AUGLAG, nlopt.LN_AUGLAG_EQ\n",
    "\n",
    "# if nlopt.AUGLAG, we have to define a subsidiary algo\n",
    "suboptimizer = nlopt.LD_SLSQP #nlopt.LN_COBYLA\n",
    "\n",
    "# define objective function\n",
    "def f(x, grad):\n",
    "    global n_iter\n",
    "    n_iter += 1\n",
    "    if grad.size > 0:\n",
    "        grad[:] = 2*x.T.dot(C)\n",
    "    return x.T.dot(C).dot(x)\n",
    "\n",
    "# define norm constraint\n",
    "def c1(x, grad):\n",
    "    if grad.size > 0:\n",
    "        grad[:] = 2*x\n",
    "    return (x.T.dot(x) - 1)\n",
    "\n",
    "# create optimizer\n",
    "n = X.shape[1] # nb of parameters to optimize, size of the vector\n",
    "opt = nlopt.opt(optimizer, n)\n",
    "print(\"Algo: %s\" % opt.get_algorithm_name())\n",
    "opt.set_max_objective(f)\n",
    "\n",
    "# if nlopt.GN_ISRES, we can define the population size\n",
    "opt.set_population(0) # by default for ISRES: pop=20*(n+1)\n",
    "\n",
    "# if nlopt.AUGLAG, set the subsidiary algo\n",
    "subopt = nlopt.opt(suboptimizer, n)\n",
    "subopt.set_lower_bounds(-1)\n",
    "subopt.set_upper_bounds(1)\n",
    "#subopt.set_ftol_rel(1e-2)\n",
    "#subopt.set_maxeval(100)\n",
    "opt.set_local_optimizer(subopt)\n",
    "\n",
    "# define bound constraints (should be between -1 and 1 because the norm should be 1)\n",
    "opt.set_lower_bounds(-1.)\n",
    "opt.set_upper_bounds(1.)\n",
    "\n",
    "# define equality constraints\n",
    "opt.add_equality_constraint(c1, 0)\n",
    "#opt.add_equality_mconstraint(constraints, tol)\n",
    "\n",
    "# define stopping criteria\n",
    "#opt.set_stopval(stopval)\n",
    "opt.set_ftol_rel(1e-8)\n",
    "#opt.set_xtol_rel(1e-4)\n",
    "opt.set_maxeval(100000) # nb of iteration\n",
    "opt.set_maxtime(10) # time in secs\n",
    "\n",
    "# define initial value\n",
    "#x0 = np.zeros((n,))\n",
    "x0 = np.array([0.1,0.1])\n",
    "\n",
    "# optimize\n",
    "x = x0\n",
    "try:\n",
    "    x = opt.optimize(x0)\n",
    "except nlopt.RoundoffLimited as e:\n",
    "    pass\n",
    "max_value = opt.last_optimum_value()\n",
    "result = opt.last_optimize_result()\n",
    "\n",
    "print(\"Max value: %f\" % max_value)\n",
    "print(\"Nb of iterations: %d\" % n_iter)\n",
    "print(\"Result: %s\" % nlopt_results[result])\n",
    "print(\"Optimized array:\")\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define PCA optimization method formally using NLopt\n",
    "def center_data(X, normalize=False, copy=True):\n",
    "    if copy:\n",
    "        X = np.copy(X)\n",
    "\n",
    "    # center the data\n",
    "    mean = X.mean(axis=0)\n",
    "    X -= mean\n",
    "    N = X.shape[0]\n",
    "    T = X.shape[1]\n",
    "\n",
    "    # normalize\n",
    "    if normalize:\n",
    "        X /= X.std(axis=0)\n",
    "    \n",
    "    return X\n",
    "\n",
    "class OrthogonalConstraint(object):\n",
    "    \n",
    "    def __init__(self, v):\n",
    "        self.v = np.copy(v)\n",
    "        \n",
    "    def constraint(self, x, grad):\n",
    "        if grad.size > 0:\n",
    "            grad[:] = self.v\n",
    "        return (x.T.dot(self.v))\n",
    "        \n",
    "\n",
    "def pca_nlopt(X, method=None, submethod=None, rough_param=0.0, normalize=False, copy=True):\n",
    "    \"\"\"\n",
    "    Compute PCA on the given data using nlopt.\n",
    "    \n",
    "    :param (str) method: it can take the following value: 'SLSQP', 'ISRES',\n",
    "                         'COBYLA', 'AUGLAG'. By default, it will be 'SLSQP'.\n",
    "    :param (str) submethod: this needs to be defined only if method is 'AUGLAG'.\n",
    "                            By default, it will be 'SLSQP'.\n",
    "    \"\"\"\n",
    "    # center the data\n",
    "    X = center_data(X, normalize=normalize, copy=copy)\n",
    "\n",
    "    # compute 'covariance' matrix and cache it\n",
    "    N = X.shape[0]\n",
    "    C = X.T.dot(X) / (N-1)\n",
    "    \n",
    "    # define useful variables\n",
    "    nlopt_results = {1: 'success', 2: 'stop_val reached', 3: 'ftol reached', 4: 'xtol reached',\n",
    "                 5: 'maxeval reached', 6: 'maxtime reached', -1: 'failure', -2: 'invalid args',\n",
    "                 -3: 'out of memory', -4: 'roundoff limited', -5: 'forced stop'}\n",
    "    n = X.shape[1] # nb of parameters to optimize, size of the vector\n",
    "    \n",
    "    # define random seed\n",
    "    nlopt.srand(125)\n",
    "\n",
    "    # define which solver (and possibly subsolver) to use\n",
    "    def get_opt(method):\n",
    "        if method == 'ISRES':\n",
    "            return nlopt.opt(nlopt.GN_ISRES, n)\n",
    "        elif method == 'COBYLA':\n",
    "            return nlopt.opt(nlopt.LN_COBYLA, n)\n",
    "        elif method == 'SLSQP':\n",
    "            return nlopt.opt(nlopt.LD_SLSQP, n)\n",
    "        elif method == 'AUGLAG':\n",
    "            return nlopt.opt(nlopt.AUGLAG, n)\n",
    "        else:\n",
    "            raise NotImplementedError(\"The given method has not been implemented\")\n",
    "\n",
    "    if method is None:\n",
    "        method = 'SLSQP'            \n",
    "    opt = get_opt(method)\n",
    "    if method == 'AUGLAG':\n",
    "        if submethod is None:\n",
    "            submethod = 'SLSQP'\n",
    "        elif submethod == 'AUGLAG':\n",
    "            raise ValueError(\"Submethod should be different from AUGLAG\")\n",
    "        subopt = get_opt(submethod)\n",
    "        subopt.set_lower_bounds(-1)\n",
    "        subopt.set_upper_bounds(1)\n",
    "        #subopt.set_ftol_rel(1e-2)\n",
    "        #subopt.set_maxeval(100)\n",
    "        opt.set_local_optimizer(subopt)\n",
    "        \n",
    "    # define objective function\n",
    "    def f(x, grad):\n",
    "        if grad.size > 0:\n",
    "            grad[:] = 2*x.T.dot(C)\n",
    "        return x.T.dot(C).dot(x)\n",
    "\n",
    "    # define norm constraint\n",
    "    def c1(x, grad):\n",
    "        if grad.size > 0:\n",
    "            grad[:] = 2*x\n",
    "        return (x.T.dot(x) - 1)\n",
    "    \n",
    "    # define objective function\n",
    "    opt.set_max_objective(f)\n",
    "        \n",
    "    # if nlopt.GN_ISRES, we can define the population size\n",
    "    opt.set_population(0) # by default for ISRES: pop=20*(n+1)\n",
    "    \n",
    "    # define bound constraints (should be between -1 and 1 because the norm should be 1)\n",
    "    opt.set_lower_bounds(-1.)\n",
    "    opt.set_upper_bounds(1.)\n",
    "\n",
    "    # define equality constraints\n",
    "    opt.add_equality_constraint(c1, 0)\n",
    "    #opt.add_equality_mconstraint(constraints, tol)\n",
    "\n",
    "    # define stopping criteria\n",
    "    #opt.set_stopval(stopval)\n",
    "    opt.set_ftol_rel(1e-8)\n",
    "    #opt.set_xtol_rel(1e-4)\n",
    "    opt.set_maxeval(100000) # nb of iteration\n",
    "    opt.set_maxtime(2) # time in secs\n",
    "\n",
    "    # define initial value\n",
    "    x0 = np.array([0.1]*n) # important that the initial value ≠ 0 for the computation of the grad!\n",
    "\n",
    "    evals, evecs, msgs = [], [], {}\n",
    "    for i in range(n):\n",
    "        if i > 0:\n",
    "            c = OrthogonalConstraint(x)\n",
    "            opt.add_equality_constraint(c.constraint, 0)\n",
    "        # optimize\n",
    "        try:\n",
    "            x = opt.optimize(x0)\n",
    "        except nlopt.RoundoffLimited as e:\n",
    "            pass\n",
    "\n",
    "        # save values\n",
    "        evecs.append(x)\n",
    "        evals.append(opt.last_optimum_value())\n",
    "        msgs[i] = nlopt_results[opt.last_optimize_result()]\n",
    "\n",
    "    return np.array(evals), np.array(evecs).T, msgs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "method = 'SLSQP' # 'SLSQP', 'COBYLA', 'ISRES', 'AUGLAG'\n",
    "submethod = None\n",
    "\n",
    "evals, P, msgs = pca_nlopt(X, method=method, submethod=submethod)\n",
    "print(msgs)\n",
    "print(evals)\n",
    "print(P)\n",
    "\n",
    "# Plot\n",
    "plt.figure(figsize=(10,4))\n",
    "plt.subplot(1,2,1)\n",
    "plt.title('PCA: eigenvalues')\n",
    "plt.bar(np.array([0.,0.1]), evals, width=0.1)\n",
    "plt.xlim(0.,1.)\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.title('PCA: data and eigenvectors')\n",
    "plt.scatter(X[:,0], X[:,1], color='b')\n",
    "plt.arrow(0, 0, P[0,0], P[1,0], length_includes_head = True, head_width = 0.15, color='b')\n",
    "plt.arrow(0, 0, P[0,1], P[1,1], length_includes_head = True, head_width = 0.15, color='g')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For comparison purpose, we obtained the following values with scipy ('SLSQP'):\n",
    "\n",
    "[ 7.93954312  0.06045688] <br>\n",
    "[[ 0.83849224 -0.54491355] <br>\n",
    " [ 0.54491353  0.83849226]]\n",
    "\n",
    "And these values using std PCA:\n",
    "\n",
    "[ 7.93954312  0.06045688] <br>\n",
    "[[-0.83849224  0.54491354] <br>\n",
    " [-0.54491354 -0.83849224]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using IPopt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Playground with IPopt\n",
    "import ipopt\n",
    "\n",
    "# define useful vars\n",
    "n = X.shape[1]\n",
    "N = X.shape[0]\n",
    "C = X.T.dot(X) / (N-1)\n",
    "\n",
    "# define initial value\n",
    "x0 = np.array([0.1]*n)\n",
    "\n",
    "# define (lower and upper) bound constraints\n",
    "lb = [-1]*n\n",
    "ub = [1]*n\n",
    "\n",
    "# define constraints\n",
    "cl = [1]\n",
    "cu = [1]\n",
    "\n",
    "class Opt(object):\n",
    "    \n",
    "    def __init__(self, verbose=True):\n",
    "        self.verbose = verbose\n",
    "        self.iter_count = 0\n",
    "\n",
    "    def objective(self, x):\n",
    "        # objective fct to minimize\n",
    "        return -x.T.dot(C).dot(x)\n",
    "    \n",
    "    def gradient(self, x):\n",
    "        # grad of the objective fct\n",
    "        return -2*x.T.dot(C)\n",
    "    \n",
    "    def constraints(self, x):\n",
    "        # norm constraint\n",
    "        return x.T.dot(x)\n",
    "    \n",
    "    def jacobian(self, x):\n",
    "        return 2*x\n",
    "    \n",
    "    #def hessian(self, x):\n",
    "    #    pass\n",
    "    \n",
    "    def intermediate(self, alg_mod, iter_count, obj_value, inf_pr, inf_du, mu, d_norm,\n",
    "                     regularization_size, alpha_du, alpha_pr, ls_trials):\n",
    "        if self.verbose:\n",
    "            print(\"Objective value at iteration #%d: %g\" % (iter_count, obj_value))\n",
    "        self.iter_count = iter_count\n",
    "\n",
    "opt = Opt(verbose=False)\n",
    "nlp = ipopt.problem(n=n, m=len(cl), problem_obj=opt, lb=lb, ub=ub, cl=cl, cu=cu)\n",
    "\n",
    "x, info = nlp.solve(x0)\n",
    "print(\"Max value: %f\" % -info['obj_val'])\n",
    "print(\"Nb of iterations: %d\" % opt.iter_count)\n",
    "print(\"Result: %s\" % info['status_msg'])\n",
    "print(\"Optimized array:\")\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define PCA optimization method formally using IPopt\n",
    "class NormConstraint(object):\n",
    "    \n",
    "    def __init__(self):\n",
    "        pass\n",
    "    \n",
    "    def constraint(self, x):\n",
    "        return x.T.dot(x)\n",
    "    \n",
    "    def jacobian(self, x):\n",
    "        return 2*x\n",
    "    \n",
    "class OrthogonalConstraint(object):\n",
    "    \n",
    "    def __init__(self, v):\n",
    "        self.v = np.copy(v)\n",
    "    \n",
    "    def constraint(self, x):\n",
    "        return x.T.dot(self.v)\n",
    "    \n",
    "    def jacobian(self, x):\n",
    "        return self.v\n",
    "\n",
    "class IPopt(object):\n",
    "\n",
    "    def __init__(self, verbose=True):\n",
    "        self.verbose = verbose\n",
    "        self.iter_count = 0\n",
    "        self.cons = []\n",
    "    \n",
    "    def add_constraint(self, c):\n",
    "        self.cons.append(c)\n",
    "\n",
    "    def objective(self, x):\n",
    "        # objective fct to minimize\n",
    "        return -x.T.dot(C).dot(x)\n",
    "    \n",
    "    def gradient(self, x):\n",
    "        # grad of the objective fct\n",
    "        return -2*x.T.dot(C)\n",
    "    \n",
    "    def constraints(self, x):\n",
    "        return np.array([c.constraint(x) for c in self.cons])\n",
    "    \n",
    "    def jacobian(self, x):\n",
    "        return np.array([c.jacobian(x) for c in self.cons])\n",
    "    \n",
    "    #def hessian(self, x):\n",
    "    #    pass\n",
    "    \n",
    "    def intermediate(self, alg_mod, iter_count, obj_value, inf_pr, inf_du, mu, d_norm,\n",
    "                     regularization_size, alpha_du, alpha_pr, ls_trials):\n",
    "        if self.verbose:\n",
    "            print(\"Objective value at iteration #%d: %g\" % (iter_count, obj_value))\n",
    "        self.iter_count = iter_count\n",
    "\n",
    "def pca_ipopt(X, rough_param=0.0, normalize=False, copy=True):\n",
    "    \"\"\"\n",
    "    Compute PCA on the given data using ipopt.\n",
    "    \"\"\"\n",
    "    # center the data\n",
    "    X = center_data(X, normalize=normalize, copy=copy)\n",
    "    \n",
    "    # define useful vars\n",
    "    n = X.shape[1]\n",
    "    N = X.shape[0]\n",
    "    C = X.T.dot(X) / (N-1)\n",
    "\n",
    "    # define initial value\n",
    "    x0 = np.array([0.1]*n)\n",
    "\n",
    "    # define (lower and upper) bound constraints\n",
    "    lb = [-1]*n\n",
    "    ub = [1]*n\n",
    "\n",
    "    # define constraints\n",
    "    cl = [1] + [0]*(n-1)\n",
    "    cu = [1] + [0]*(n-1)\n",
    "\n",
    "    evals, evecs, msgs = [], [], {}\n",
    "    opt = IPopt(verbose=False)\n",
    "    opt.add_constraint(NormConstraint())\n",
    "    for i in range(n):\n",
    "        if i > 0:\n",
    "            opt.add_constraint(OrthogonalConstraint(x))\n",
    "        \n",
    "        i1 = i+1\n",
    "        nlp = ipopt.problem(n=n, m=len(cl[:i1]), problem_obj=opt,\n",
    "                            lb=lb, ub=ub, cl=cl[:i1], cu=cu[:i1])\n",
    "        \n",
    "        # solve problem\n",
    "        x, info = nlp.solve(x0)\n",
    "        \n",
    "        evecs.append(x)\n",
    "        evals.append(-info['obj_val'])\n",
    "        msgs[i] = info['status_msg']\n",
    "        \n",
    "    return np.array(evals), np.array(evecs).T, msgs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evals, P, msgs = pca_ipopt(X)\n",
    "print(msgs)\n",
    "print(evals)\n",
    "print(P)\n",
    "\n",
    "# Plot\n",
    "plt.figure(figsize=(10,4))\n",
    "plt.subplot(1,2,1)\n",
    "plt.title('PCA: eigenvalues')\n",
    "plt.bar(np.array([0.,0.1]), evals, width=0.1)\n",
    "plt.xlim(0.,1.)\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.title('PCA: data and eigenvectors')\n",
    "plt.scatter(X[:,0], X[:,1], color='b')\n",
    "plt.arrow(0, 0, P[0,0], P[1,0], length_includes_head = True, head_width = 0.15, color='b')\n",
    "plt.arrow(0, 0, P[0,1], P[1,1], length_includes_head = True, head_width = 0.15, color='g')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}