{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Faces dataset decompositions\n\nThis example applies to `olivetti_faces_dataset` different unsupervised\nmatrix decomposition (dimension reduction) methods from the module\n:mod:`sklearn.decomposition` (see the documentation chapter\n`decompositions`).\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: The scikit-learn developers\n# SPDX-License-Identifier: BSD-3-Clause"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Dataset preparation\n\nLoading and preprocessing the Olivetti faces dataset.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import logging\n\nimport matplotlib.pyplot as plt\nfrom numpy.random import RandomState\n\nfrom sklearn import cluster, decomposition\nfrom sklearn.datasets import fetch_olivetti_faces\n\nrng = RandomState(0)\n\n# Display progress logs on stdout\nlogging.basicConfig(level=logging.INFO, format=\"%(asctime)s %(levelname)s %(message)s\")\n\nfaces, _ = fetch_olivetti_faces(return_X_y=True, shuffle=True, random_state=rng)\nn_samples, n_features = faces.shape\n\n# Global centering (focus on one feature, centering all samples)\nfaces_centered = faces - faces.mean(axis=0)\n\n# Local centering (focus on one sample, centering all features)\nfaces_centered -= faces_centered.mean(axis=1).reshape(n_samples, -1)\n\nprint(\"Dataset consists of %d faces\" % n_samples)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define a base function to plot the gallery of faces.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n_row, n_col = 2, 3\nn_components = n_row * n_col\nimage_shape = (64, 64)\n\n\ndef plot_gallery(title, images, n_col=n_col, n_row=n_row, cmap=plt.cm.gray):\n    fig, axs = plt.subplots(\n        nrows=n_row,\n        ncols=n_col,\n        figsize=(2.0 * n_col, 2.3 * n_row),\n        facecolor=\"white\",\n        constrained_layout=True,\n    )\n    fig.set_constrained_layout_pads(w_pad=0.01, h_pad=0.02, hspace=0, wspace=0)\n    fig.set_edgecolor(\"black\")\n    fig.suptitle(title, size=16)\n    for ax, vec in zip(axs.flat, images):\n        vmax = max(vec.max(), -vec.min())\n        im = ax.imshow(\n            vec.reshape(image_shape),\n            cmap=cmap,\n            interpolation=\"nearest\",\n            vmin=-vmax,\n            vmax=vmax,\n        )\n        ax.axis(\"off\")\n\n    fig.colorbar(im, ax=axs, orientation=\"horizontal\", shrink=0.99, aspect=40, pad=0.01)\n    plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's take a look at our data. Gray color indicates negative values,\nwhite indicates positive values.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plot_gallery(\"Faces from dataset\", faces_centered[:n_components])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Decomposition\n\nInitialise different estimators for decomposition and fit each\nof them on all images and plot some results. Each estimator extracts\n6 components as vectors $h \\in \\mathbb{R}^{4096}$.\nWe just displayed these vectors in human-friendly visualisation as 64x64 pixel images.\n\nRead more in the `User Guide <decompositions>`.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Eigenfaces - PCA using randomized SVD\nLinear dimensionality reduction using Singular Value Decomposition (SVD) of the data\nto project it to a lower dimensional space.\n\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>The Eigenfaces estimator, via the :py:mod:`sklearn.decomposition.PCA`,\n    also provides a scalar `noise_variance_` (the mean of pixelwise variance)\n    that cannot be displayed as an image.</p></div>\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pca_estimator = decomposition.PCA(\n    n_components=n_components, svd_solver=\"randomized\", whiten=True\n)\npca_estimator.fit(faces_centered)\nplot_gallery(\n    \"Eigenfaces - PCA using randomized SVD\", pca_estimator.components_[:n_components]\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Non-negative components - NMF\n\nEstimate non-negative original data as production of two non-negative matrices.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nmf_estimator = decomposition.NMF(n_components=n_components, tol=5e-3)\nnmf_estimator.fit(faces)  # original non- negative dataset\nplot_gallery(\"Non-negative components - NMF\", nmf_estimator.components_[:n_components])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Independent components - FastICA\nIndependent component analysis separates a multivariate vectors into additive\nsubcomponents that are maximally independent.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ica_estimator = decomposition.FastICA(\n    n_components=n_components, max_iter=400, whiten=\"arbitrary-variance\", tol=15e-5\n)\nica_estimator.fit(faces_centered)\nplot_gallery(\n    \"Independent components - FastICA\", ica_estimator.components_[:n_components]\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Sparse components - MiniBatchSparsePCA\n\nMini-batch sparse PCA (:class:`~sklearn.decomposition.MiniBatchSparsePCA`)\nextracts the set of sparse components that best reconstruct the data. This\nvariant is faster but less accurate than the similar\n:class:`~sklearn.decomposition.SparsePCA`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "batch_pca_estimator = decomposition.MiniBatchSparsePCA(\n    n_components=n_components, alpha=0.1, max_iter=100, batch_size=3, random_state=rng\n)\nbatch_pca_estimator.fit(faces_centered)\nplot_gallery(\n    \"Sparse components - MiniBatchSparsePCA\",\n    batch_pca_estimator.components_[:n_components],\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Dictionary learning\n\nBy default, :class:`~sklearn.decomposition.MiniBatchDictionaryLearning`\ndivides the data into mini-batches and optimizes in an online manner by\ncycling over the mini-batches for the specified number of iterations.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "batch_dict_estimator = decomposition.MiniBatchDictionaryLearning(\n    n_components=n_components, alpha=0.1, max_iter=50, batch_size=3, random_state=rng\n)\nbatch_dict_estimator.fit(faces_centered)\nplot_gallery(\"Dictionary learning\", batch_dict_estimator.components_[:n_components])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Cluster centers - MiniBatchKMeans\n\n:class:`sklearn.cluster.MiniBatchKMeans` is computationally efficient and\nimplements on-line learning with a\n:meth:`~sklearn.cluster.MiniBatchKMeans.partial_fit` method. That is\nwhy it could be beneficial to enhance some time-consuming algorithms with\n:class:`~sklearn.cluster.MiniBatchKMeans`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "kmeans_estimator = cluster.MiniBatchKMeans(\n    n_clusters=n_components,\n    tol=1e-3,\n    batch_size=20,\n    max_iter=50,\n    random_state=rng,\n)\nkmeans_estimator.fit(faces_centered)\nplot_gallery(\n    \"Cluster centers - MiniBatchKMeans\",\n    kmeans_estimator.cluster_centers_[:n_components],\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Factor Analysis components - FA\n\n:class:`~sklearn.decomposition.FactorAnalysis` is similar to\n:class:`~sklearn.decomposition.PCA` but has the advantage of modelling the\nvariance in every direction of the input space independently (heteroscedastic\nnoise). Read more in the `User Guide <FA>`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fa_estimator = decomposition.FactorAnalysis(n_components=n_components, max_iter=20)\nfa_estimator.fit(faces_centered)\nplot_gallery(\"Factor Analysis (FA)\", fa_estimator.components_[:n_components])\n\n# --- Pixelwise variance\nplt.figure(figsize=(3.2, 3.6), facecolor=\"white\", tight_layout=True)\nvec = fa_estimator.noise_variance_\nvmax = max(vec.max(), -vec.min())\nplt.imshow(\n    vec.reshape(image_shape),\n    cmap=plt.cm.gray,\n    interpolation=\"nearest\",\n    vmin=-vmax,\n    vmax=vmax,\n)\nplt.axis(\"off\")\nplt.title(\"Pixelwise variance from \\n Factor Analysis (FA)\", size=16, wrap=True)\nplt.colorbar(orientation=\"horizontal\", shrink=0.8, pad=0.03)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Decomposition: Dictionary learning\n\nIn the further section, let's consider `DictionaryLearning` more precisely.\nDictionary learning is a problem that amounts to finding a sparse representation\nof the input data as a combination of simple elements. These simple elements form\na dictionary. It is possible to constrain the dictionary and/or coding coefficients\nto be positive to match constraints that may be present in the data.\n\n:class:`~sklearn.decomposition.MiniBatchDictionaryLearning` implements a\nfaster, but less accurate version of the dictionary learning algorithm that\nis better suited for large datasets. Read more in the `User Guide\n<MiniBatchDictionaryLearning>`.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the same samples from our dataset but with another colormap.\nRed indicates negative values, blue indicates positive values,\nand white represents zeros.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plot_gallery(\"Faces from dataset\", faces_centered[:n_components], cmap=plt.cm.RdBu)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Similar to the previous examples, we change parameters and train\n:class:`~sklearn.decomposition.MiniBatchDictionaryLearning` estimator on all\nimages. Generally, the dictionary learning and sparse encoding decompose\ninput data into the dictionary and the coding coefficients matrices. $X\n\\approx UV$, where $X = [x_1, . . . , x_n]$, $X \\in\n\\mathbb{R}^{m\u00d7n}$, dictionary $U \\in \\mathbb{R}^{m\u00d7k}$, coding\ncoefficients $V \\in \\mathbb{R}^{k\u00d7n}$.\n\nAlso below are the results when the dictionary and coding\ncoefficients are positively constrained.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Dictionary learning - positive dictionary\n\nIn the following section we enforce positivity when finding the dictionary.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dict_pos_dict_estimator = decomposition.MiniBatchDictionaryLearning(\n    n_components=n_components,\n    alpha=0.1,\n    max_iter=50,\n    batch_size=3,\n    random_state=rng,\n    positive_dict=True,\n)\ndict_pos_dict_estimator.fit(faces_centered)\nplot_gallery(\n    \"Dictionary learning - positive dictionary\",\n    dict_pos_dict_estimator.components_[:n_components],\n    cmap=plt.cm.RdBu,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Dictionary learning - positive code\n\nBelow we constrain the coding coefficients as a positive matrix.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dict_pos_code_estimator = decomposition.MiniBatchDictionaryLearning(\n    n_components=n_components,\n    alpha=0.1,\n    max_iter=50,\n    batch_size=3,\n    fit_algorithm=\"cd\",\n    random_state=rng,\n    positive_code=True,\n)\ndict_pos_code_estimator.fit(faces_centered)\nplot_gallery(\n    \"Dictionary learning - positive code\",\n    dict_pos_code_estimator.components_[:n_components],\n    cmap=plt.cm.RdBu,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Dictionary learning - positive dictionary & code\n\nAlso below are the results if the dictionary values and coding\ncoefficients are positively constrained.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dict_pos_estimator = decomposition.MiniBatchDictionaryLearning(\n    n_components=n_components,\n    alpha=0.1,\n    max_iter=50,\n    batch_size=3,\n    fit_algorithm=\"cd\",\n    random_state=rng,\n    positive_dict=True,\n    positive_code=True,\n)\ndict_pos_estimator.fit(faces_centered)\nplot_gallery(\n    \"Dictionary learning - positive dictionary & code\",\n    dict_pos_estimator.components_[:n_components],\n    cmap=plt.cm.RdBu,\n)"
      ]
    }
  ],
  "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.9.21"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}