{ "cells": [ { "cell_type": "raw", "metadata": {}, "source": [ "---\n", "title: Biostat/Biomath M257 Homework 2\n", "subtitle: 'Due Apr 28 @ 11:59PM'\n", "author: Student Name and UID\n", "date: today\n", "format:\n", " html:\n", " theme: cosmo\n", " embed-resources: true\n", " number-sections: true\n", " toc: true\n", " toc-depth: 4\n", " toc-location: left\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "System information (for reproducibility):" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Julia Version 1.8.5\n", "Commit 17cfb8e65ea (2023-01-08 06:45 UTC)\n", "Platform Info:\n", " OS: macOS (arm64-apple-darwin21.5.0)\n", " CPU: 12 × Apple M2 Max\n", " WORD_SIZE: 64\n", " LIBM: libopenlibm\n", " LLVM: libLLVM-13.0.1 (ORCJIT, apple-m1)\n", " Threads: 8 on 8 virtual cores\n", "Environment:\n", " JULIA_NUM_THREADS = 8\n", " JULIA_EDITOR = code\n" ] } ], "source": [ "versioninfo()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Load packages:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "tags": [] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[32m\u001b[1m Activating\u001b[22m\u001b[39m project at `~/Documents/github.com/ucla-biostat-257/2023spring/hw/hw2`\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\u001b[32m\u001b[1mStatus\u001b[22m\u001b[39m `~/Documents/github.com/ucla-biostat-257/2023spring/hw/hw2/Project.toml`\n", " \u001b[90m [6e4b80f9] \u001b[39mBenchmarkTools v1.3.2\n", " \u001b[90m [7073ff75] \u001b[39mIJulia v1.24.0\n", " \u001b[90m [916415d5] \u001b[39mImages v0.25.2\n", " \u001b[90m [bdcacae8] \u001b[39mLoopVectorization v0.12.157\n", " \u001b[90m [8bb1440f] \u001b[39mDelimitedFiles\n", " \u001b[90m [37e2e46d] \u001b[39mLinearAlgebra\n", " \u001b[90m [9abbd945] \u001b[39mProfile\n" ] } ], "source": [ "using Pkg\n", "\n", "Pkg.activate(pwd())\n", "Pkg.instantiate()\n", "Pkg.status()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "tags": [] }, "outputs": [], "source": [ "# load libraries\n", "using BenchmarkTools, DelimitedFiles, Images, LinearAlgebra, LoopVectorization\n", "using Profile, Random" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Q1. Nonnegative Matrix Factorization\n", "\n", "Nonnegative matrix factorization (NNMF) was introduced by [Lee and Seung (1999)](https://www.nature.com/articles/44565) as an alternative to principal components and vector quantization with applications in data compression, clustering, and deconvolution. In this homework we consider algorithms for fitting NNMF and (optionally) high performance computing using graphical processing units (GPUs).\n", "\n", "\n", "\n", "In mathematical terms, one approximates a data matrix $\\mathbf{X} \\in \\mathbb{R}^{m \\times n}$ with nonnegative entries $x_{ij}$ by a product of two low-rank matrices $\\mathbf{V} \\in \\mathbb{R}^{m \\times r}$ and $\\mathbf{W} \\in \\mathbb{R}^{r \\times n}$ with nonnegative entries $v_{ik}$ and $w_{kj}$. Consider minimization of the squared Frobenius norm\n", "$$\n", "L(\\mathbf{V}, \\mathbf{W}) = \\|\\mathbf{X} - \\mathbf{V} \\mathbf{W}\\|_{\\text{F}}^2 = \\sum_i \\sum_j \\left(x_{ij} - \\sum_k v_{ik} w_{kj} \\right)^2, \\quad v_{ik} \\ge 0, w_{kj} \\ge 0,\n", "$$\n", "which should lead to a good factorization. Lee and Seung suggest an iterative algorithm with multiplicative updates\n", "$$\n", "v_{ik}^{(t+1)} = v_{ik}^{(t)} \\frac{\\sum_j x_{ij} w_{kj}^{(t)}}{\\sum_j b_{ij}^{(t)} w_{kj}^{(t)}}, \\quad \\text{where } b_{ij}^{(t)} = \\sum_k v_{ik}^{(t)} w_{kj}^{(t)},\n", "$$\n", "$$\n", "w_{kj}^{(t+1)} = w_{kj}^{(t)} \\frac{\\sum_i x_{ij} v_{ik}^{(t+1)}}{\\sum_i b_{ij}^{(t+1/2)} v_{ik}^{(t+1)}}, \\quad \\text{where } b_{ij}^{(t+1/2)} = \\sum_k v_{ik}^{(t+1)} w_{kj}^{(t)}\n", "$$\n", "that will drive the objective $L^{(t)} = L(\\mathbf{V}^{(t)}, \\mathbf{W}^{(t)})$ downhill. Superscript $t$ indicates the iteration number. In following questions, efficiency (both speed and memory) will be the most important criterion when grading this problem." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q1.1 Develop code\n", "\n", "Implement the algorithm with arguments: $\\mathbf{X}$ (data, each row is a vectorized image), rank $r$, convergence tolerance, and optional starting point." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "nnmf (generic function with 1 method)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function nnmf(\n", " # positional arguments\n", " X :: AbstractMatrix{T}, \n", " r :: Integer;\n", " # kw arguments\n", " maxiter :: Integer = 1000, \n", " tolfun :: Number = 1e-4,\n", " V :: AbstractMatrix{T} = Random.rand!(similar(X, size(X, 1), r)),\n", " W :: AbstractMatrix{T} = Random.rand!(similar(X, r, size(X, 2))),\n", " ) where T <: AbstractFloat\n", " # TODO: implementation\n", " # Output\n", " V, W, obj, niter\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q1.2 Data\n", "\n", "Database 1 from the [MIT Center for Biological and Computational Learning (CBCL)](http://cbcl.mit.edu/software-datasets/FaceData2.html) reduces to a matrix $\\mathbf{X}$ containing $m = 2,429$ gray-scale face images with $n = 19 \\times 19 = 361$ pixels per face. Each image (row) is scaled to have mean and standard deviation 0.25. \n", "\n", "Read in the [`nnmf-2429-by-361-face.txt`](https://raw.githubusercontent.com/ucla-biostat-257/2023spring/master/hw/hw2/nnmf-2429-by-361-face.txt) file, e.g., using [`readdlm`](https://docs.julialang.org/en/v1/stdlib/DelimitedFiles/#Delimited-Files) function, and display a couple of sample images, e.g., using the [Images.jl](https://juliaimages.org/stable/) package." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHIAAAByCAAAAACqttqhAAAABGdBTUEAALGPC/xhBQAAAAFzUkdCAK7OHOkAAAAgY0hSTQAAeiYAAICEAAD6AAAAgOgAAHUwAADqYAAAOpgAABdwnLpRPAAAAydJREFUaAW9wTFr3AUcx+EPyfc40oZrSogNh1JCQzlaQiQURdxcRHDwPbg7ODv5IsQX4OYglIK4SKGUaildQksJxnA2hMbjb45wzXmcCfKZ/lA6/54nN9Et9D4aoC6aojV0Da3SOkb30Sv0GDXoJtpFMzREoVwoF8qFcqFcDtAcnaInaA816CO0gm6iLXSIvqa1gsZolzd1USgXyoVyoVwolwV0iA7QBjqgdQ99hu6hL9B3tNZRQ2sZTVAXzVAoF8qFcqFcKJcZb/oSfYMGaIaeo2X0M/oX7aCP0bfoKprQmqEFFMqFcqFcKBfKZQ2NaH2PtlEH7aPL6IzWFfQANWiC+ugETWl1USgXyoVyoVwolw7qoVM0pTVBPXQD7aMFtERriO6iFbSMztCEVigXyoVyoVwol0/RI7SN+miGrqC/UR/9hUZoiHpoFe2gPmrQErpPK5QL5UK5UC6Uyw00RAN0HfXRHN1FY3SOOmgVnaItWhtoDY3QOlpFoVwoF8qFcqFcRuhzdIT6aBFdoD7aRRdogIa0eugW6qAtdILO0WsUyoVyoVwoF8rlFB2gd1EH9dAGOkO/oAG6jX5FHbSKPkTnaIymqEMrlAvlQrlQLpTLCK2iI3SBVtAAvYOOUR+toG30AzpCW+h3NEcNGqMjFMqFcqFcKBfK5QyN0BJq0Bw9RItoAT1AX6E9NEdj9BPaRGPUoAUUFMqFcqFcKBfKpYeCZqiPNtAYvUTvoSE6QMu0XqIf0Q66gy6hGa1QLpQL5UK5UC6HqIsW0RQ9Ro/QM3QJbaOr6BxtoT30FL1AD9EdNEHHKJQL5UK5UC6UyxCtoT56ie6jMdpEAzRAi+gE9dAMnaFDdI56aAMNUSgXyoVyoVwolzFq0FV0gP5DE7SLdtES2kR/oClv9xpdQ6eoQaFcKBfKhXKhXC6j52gdraGnvN0UvUBz3vQJmqITdBs9oBXKhXKhXCgXyuU1mqMRuoE2UBddQmfoOlpAY9RBe2gN/YNm6BX6k1YoF8qFcqFcKJc5rSfoGdpEH6ARalCDJqiLGnQH7aN1dIJ+Q/u0QrlQLpQL5UK5/wFd2bFZ/j51mAAAAABJRU5ErkJggg==", "text/html": [ "" ], "text/plain": [ "19×19 reinterpret(reshape, Gray{Float64}, ::Matrix{Float64}) with eltype Gray{Float64}:\n", " Gray{Float64}(0.14815) … Gray{Float64}(0.0)\n", " Gray{Float64}(0.018294) Gray{Float64}(0.027569)\n", " Gray{Float64}(0.027569) Gray{Float64}(0.0)\n", " Gray{Float64}(0.0) Gray{Float64}(0.0090181)\n", " Gray{Float64}(0.083222) Gray{Float64}(0.036844)\n", " Gray{Float64}(0.10177) … Gray{Float64}(0.064671)\n", " Gray{Float64}(0.38004) Gray{Float64}(0.13887)\n", " Gray{Float64}(0.51917) Gray{Float64}(0.30583)\n", " Gray{Float64}(0.43569) Gray{Float64}(0.37076)\n", " Gray{Float64}(0.38004) Gray{Float64}(0.37076)\n", " Gray{Float64}(0.29656) … Gray{Float64}(0.26873)\n", " Gray{Float64}(0.25946) Gray{Float64}(0.28728)\n", " Gray{Float64}(0.31511) Gray{Float64}(0.28728)\n", " Gray{Float64}(0.2038) Gray{Float64}(0.19453)\n", " Gray{Float64}(0.083222) Gray{Float64}(0.073946)\n", " Gray{Float64}(0.018294) … Gray{Float64}(0.0)\n", " Gray{Float64}(0.064671) Gray{Float64}(0.0)\n", " Gray{Float64}(0.018294) Gray{Float64}(0.0)\n", " Gray{Float64}(0.0) Gray{Float64}(0.0)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X = readdlm(\"nnmf-2429-by-361-face.txt\")\n", "colorview(Gray, reshape(X[1, :], 19, 19))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHIAAAByCAAAAACqttqhAAAABGdBTUEAALGPC/xhBQAAAAFzUkdCAK7OHOkAAAAgY0hSTQAAeiYAAICEAAD6AAAAgOgAAHUwAADqYAAAOpgAABdwnLpRPAAAAvNJREFUaAW9wUFrmwUcx/Fvn/6aNIkxTTWmdq4yEYXJDp58FR59mV49+hbEi4gIEyvMbY40jWm7p2kq31weKMPj//MJO4/RB2iCBmiKJmiIxqih00eHqIf20RY16BZdoVAulAvlQrlQLsco6A36g84UfYLO0Bb1UQ816BY1aI9Og+7QAQrlQrlQLpQL5dJH52jDQxfoOZ2n6Gs0QdfoN/QTukJ7aB/N0TMUyoVyoVwoF8olvNseuuehX1CLnqDn6E/U0hmjE7RAP6JQLpQL5UK5UC4T9Bpt6Nzzf95Dx2gP/c5Dl+gGTdEchXKhXCgXyoVyGaDw0BG6QHvons4MPUU/o2O0QGN0iQ7RWzRBoVwoF8qFcqFchugMLdES9dGn6CW6QT20jzZojM5QHw3RF+iEzgiFcqFcKBfKhXIZoS36GG3RFt2gBi3QAbpDN6iPvkKn6BodoCmaoxEK5UK5UC6UC+XSoBYN0DGdC7RGLRqgPrpF92iE5miAggaoh7YolAvlQrlQLpRLg5Z0GhS0obOlM0YNOkBLNEFDdIBCZ4NCuVAulAvlQrk8Q7+iBeqhoBVaojfoBVqhIdpHT9Bb9CUao5bOCoVyoVwoF8qFcmHnETpH12gfXaJr1NJ5iYZogw7RMfobtWiINmiJQrlQLpQL5UK5fINeoxW6Qj10TecINWiMZugUjdAUrdE/aIpatEChXCgXyoVyoVy+RT+gY/QvukIt2kdHaI5mqEEfoVu0QA0KukFLtEahXCgXyoVyoVzYadD7qEELdEdnhdZoixq0RB/SuUUbOmvUolAulAvlQrlQLux8jl6hGXqFLlHLQy94txn6DM1QUIuCQrlQLpQL5UK5sDNDLZqgARqjc7RBQXM0orNGp+gROkAbFDRCoVwoF8qFcqFc2PkOfY/6aIxWaI5WaIQmaIRmaIbGqEFLtEINCgrlQrlQLpQL5cJOHwVt0Qy9Rj3UR7foL3SCHqMpCrpCW7RBQdcolAvlQrlQLpT7DxQjibc/vHJgAAAAAElFTkSuQmCC", "text/html": [ "" ], "text/plain": [ "19×19 reinterpret(reshape, Gray{Float64}, ::Matrix{Float64}) with eltype Gray{Float64}:\n", " Gray{Float64}(0.0) Gray{Float64}(0.1135) … Gray{Float64}(0.54206)\n", " Gray{Float64}(0.069167) Gray{Float64}(0.12828) Gray{Float64}(0.54206)\n", " Gray{Float64}(0.098722) Gray{Float64}(0.017445) Gray{Float64}(0.28344)\n", " Gray{Float64}(0.1135) Gray{Float64}(0.017445) Gray{Float64}(0.024834)\n", " Gray{Float64}(0.17261) Gray{Float64}(0.032222) Gray{Float64}(0.1135)\n", " Gray{Float64}(0.20956) Gray{Float64}(0.047) … Gray{Float64}(0.15044)\n", " Gray{Float64}(0.2465) Gray{Float64}(0.16522) Gray{Float64}(0.34256)\n", " Gray{Float64}(0.29083) Gray{Float64}(0.23911) Gray{Float64}(0.49772)\n", " Gray{Float64}(0.29822) Gray{Float64}(0.21694) Gray{Float64}(0.54944)\n", " Gray{Float64}(0.30561) Gray{Float64}(0.15783) Gray{Float64}(0.53467)\n", " Gray{Float64}(0.48294) Gray{Float64}(0.15044) … Gray{Float64}(0.48294)\n", " Gray{Float64}(0.48294) Gray{Float64}(0.25389) Gray{Float64}(0.42383)\n", " Gray{Float64}(0.69722) Gray{Float64}(0.40167) Gray{Float64}(0.35733)\n", " Gray{Float64}(1.0) Gray{Float64}(0.40906) Gray{Float64}(0.313)\n", " Gray{Float64}(1.0) Gray{Float64}(0.41644) Gray{Float64}(0.29083)\n", " Gray{Float64}(1.0) Gray{Float64}(0.55683) … Gray{Float64}(0.30561)\n", " Gray{Float64}(1.0) Gray{Float64}(0.63811) Gray{Float64}(0.34994)\n", " Gray{Float64}(1.0) Gray{Float64}(0.97061) Gray{Float64}(0.36472)\n", " Gray{Float64}(1.0) Gray{Float64}(1.0) Gray{Float64}(0.32778)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "colorview(Gray, reshape(X[10, :], 19, 19))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q1.3 Correctness and efficiency\n", "\n", "Report the run times, using `@btime`, of your function for fitting NNMF on the MIT CBCL face data set at ranks $r=10, 20, 30, 40, 50$. For ease of comparison (and grading), please start your algorithm with the provided $\\mathbf{V}^{(0)}$ (first $r$ columns of [`V0.txt`](https://raw.githubusercontent.com/ucla-biostat-257/2023spring/master/hw/hw2/V0.txt)) and $\\mathbf{W}^{(0)}$ (first $r$ rows of [`W0.txt`](https://raw.githubusercontent.com/ucla-biostat-257/2023spring/master/hw/hw2/W0.txt)) and stopping criterion\n", "$$\n", "\\frac{|L^{(t+1)} - L^{(t)}|}{|L^{(t)}| + 1} \\le 10^{-4}.\n", "$$\n", "\n", "**Hint**: When I run the following code using my own implementation of `nnmf`\n", "```julia\n", "# provided start point\n", "V0full = readdlm(\"V0.txt\", ' ', Float64)\n", "W0full = readdlm(\"W0.txt\", ' ', Float64);\n", "\n", "# benchmarking\n", "for r in [10, 20, 30, 40, 50]\n", " println(\"r=$r\")\n", " V0 = V0full[:, 1:r]\n", " W0 = W0full[1:r, :]\n", " _, _, obj, niter = nnmf(X, r, V = V0, W = W0)\n", " @btime nnmf($X, $r, V = $V0, W = $W0) setup=(\n", " copyto!(V0, V0full[:, 1:r]), \n", " copyto!(W0, W0full[1:r, :])\n", " )\n", " println(\"obj=$obj, niter=$niter\")\n", "end\n", "```\n", "the output is\n", "```\n", "r=10\n", " 162.662 ms (9 allocations: 437.19 KiB)\n", "obj=11730.866905748058, niter=239\n", "r=20\n", " 234.293 ms (9 allocations: 875.44 KiB)\n", "obj=8497.605595863002, niter=394\n", "r=30\n", " 259.524 ms (9 allocations: 1.28 MiB)\n", "obj=6621.94596847528, niter=482\n", "r=40\n", " 289.918 ms (9 allocations: 1.72 MiB)\n", "obj=5256.866299829562, niter=581\n", "r=50\n", " 397.511 ms (10 allocations: 2.15 MiB)\n", "obj=4430.362097310877, niter=698\n", "```\n", "Due to machine differences, your run times can be different from mine but certainly can not be order of magnitude longer. Your memory allocation should be less or equal to mine." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q1.4 Non-uniqueness\n", "\n", "Choose an $r \\in \\{10, 20, 30, 40, 50\\}$ and start your algorithm from a different $\\mathbf{V}^{(0)}$ and $\\mathbf{W}^{(0)}$. Do you obtain the same objective value and $(\\mathbf{V}, \\mathbf{W})$? Explain what you find." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q1.5 Fixed point\n", "\n", "For the same $r$, start your algorithm from $v_{ik}^{(0)} = w_{kj}^{(0)} = 1$ for all $i,j,k$. Do you obtain the same objective value and $(\\mathbf{V}, \\mathbf{W})$? Explain what you find." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q1.6 Interpret NNMF result\n", "\n", "Plot the basis images (rows of $\\mathbf{W}$) at rank $r=50$. What do you find?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q1.7 GPU (optional)\n", "\n", "Investigate the GPU capabilities of Julia. Report the speed gain of your GPU code over CPU code at ranks $r=10, 20, 30, 40, 50$. Make sure to use the same starting point as in Q1.3." ] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "Julia (8 threads) 1.8.5", "language": "julia", "name": "julia-_8-threads_-1.8" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.8.5" }, "toc": { "colors": { "hover_highlight": "#DAA520", "running_highlight": "#FF0000", "selected_highlight": "#FFD700" }, "moveMenuLeft": true, "nav_menu": { "height": "87px", "width": "252px" }, "navigate_menu": true, "number_sections": true, "sideBar": true, "skip_h1_title": true, "threshold": 4, "toc_cell": false, "toc_section_display": "block", "toc_window_display": false, "widenNotebook": false } }, "nbformat": 4, "nbformat_minor": 4 }