{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Biostat 257 Homework 2\n", "\n", "**Due Apr 29 @ 11:59PM**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "versioninfo()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# load libraries\n", "using BenchmarkTools, DelimitedFiles, Images, LinearAlgebra, 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 analog of principal components and vector quantization with applications in data compression and clustering. 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", "\tL(\\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 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": null, "metadata": {}, "outputs": [], "source": [ "function nnmf(\n", " X :: Matrix{T}, \n", " r :: Integer;\n", " maxiter :: Integer = 1000, \n", " tolfun :: Number = 1e-4,\n", " V :: Matrix{T} = rand(T, size(X, 1), r),\n", " W :: Matrix{T} = rand(T, r, size(X, 2))\n", " ) where T <: AbstractFloat\n", " # 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) 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/2022spring/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 sample images, e.g., using the [Images.jl](https://juliaimages.org/stable/) package." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X = readdlm(\"nnmf-2429-by-361-face.txt\")\n", "colorview(Gray, reshape(X[1, :], 19, 19))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "colorview(Gray, reshape(X[5, :], 19, 19))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q1.3 Correctness and efficiency\n", "\n", "Report the run times, using `@time`, 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/2022spring/master/hw/hw2/V0.txt)) and $\\mathbf{W}^{(0)}$ (first $r$ rows of [`W0.txt`](https://raw.githubusercontent.com/ucla-biostat-257/2022spring/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", "for r in [10, 20, 30, 40, 50]\n", " println(\"r=$r\")\n", " V0 = V0full[:, 1:r]\n", " W0 = W0full[1:r, :]\n", " @time V, W, obj, niter = nnmf(X, r; V = V0, W = W0)\n", " println(\"obj=$obj, niter=$niter\")\n", "end\n", "```\n", "the output is\n", "```\n", "r=10\n", " 1.047598 seconds (20 allocations: 6.904 MiB)\n", "obj=11730.38800985483, niter=239\n", "r=20\n", " 1.913147 seconds (20 allocations: 7.120 MiB)\n", "obj=8497.222317850326, niter=394\n", "r=30\n", " 2.434662 seconds (20 allocations: 7.336 MiB)\n", "obj=6621.627345486279, niter=482\n", "r=40\n", " 3.424469 seconds (22 allocations: 7.554 MiB)\n", "obj=5256.663870563529, niter=581\n", "r=50\n", " 4.480342 seconds (23 allocations: 7.774 MiB)\n", "obj=4430.201581697291, niter=698\n", "```\n", "Since my laptop is about 6-7 years old, I expect to see your run time shorter than mine. 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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Q2. Estimating Kinship Matrix\n", "\n", "Consider the numerical task of estimating an $n \\times n$ kinship matrix $\\Phi$ from an $n \\times m$ genotype matrix $\\mathbf{G}$. Here $n$ is the number of individuals and $m$ is the number of genetic markers. [Lange et al](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6763373/) derived a method of moment estimator of form\n", "$$\n", " \\widehat \\Phi_{ij} = \\frac{e_{ij} - \\sum_{k=1}^m [p_k^2 + (1 - p_k)^2]}{m - \\sum_{k=1}^m [p_k^2 + (1 - p_k)^2]}, \\quad 1 \\le i, j \\le n,\n", "$$\n", "where\n", "$$\n", "\\begin{eqnarray*}\n", " e_{ij} &=& \\frac{1}{4} \\sum_{k=1}^m [g_{ik} g_{jk} + (2 - g_{ik})(2 - g_{jk})] \\\\\n", " p_k &=& \\frac {1}{2n} \\sum_{i=1}^n g_{ik}.\n", "\\end{eqnarray*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q2.1 Develop code\n", "\n", "Write a function that takes a matrix `G` as input and outputs the method of moment estimator. \n", "Make your function as efficient (both speed and memory) as possible. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "function kinship(G::Matrix{T}) where T <: AbstractFloat\n", " n, m = size(G)\n", " # TODO: your code here\n", " Φ = zeros(n, m)\n", " # output\n", " Φ\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q2.2 Correctness\n", "\n", "First let's make sure our function yields correct answer. Run your function on a fake genotype matrix" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# generate a fake genotype matrix with entries {0, 1, 2}\n", "Random.seed!(257)\n", "G = rand(0.0:2.0, 1000, 10000)\n", "Φ = kinship(G)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare the upper left $5 \\times 5$ block to what I got using my implementation\n", "\n", "```julia\n", "Φ[1:5, 1:5]\n", "```\n", "\n", "```\n", "5×5 Matrix{Float64}:\n", " 0.673584 -0.000762864 -0.00266412 0.00343992 0.00293959\n", " -0.000762864 0.665178 -0.0101691 -0.0110697 0.00223912\n", " -0.00266412 -0.0101691 0.665078 0.0102444 0.00253932\n", " 0.00343992 -0.0110697 0.0102444 0.66768 -0.0083679\n", " 0.00293959 0.00223912 0.00253932 -0.0083679 0.663777\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q2.3 Efficiency\n", "\n", "In a typical genetic data set, $n$ is at order of $10^3 \\sim 10^6$ and $m$ is at order of $10^6 \\sim 10^7$. Benchmark your function using the smaller data set $G$ generated in Q2.2. Efficiency (both speed and memory) will be the most important criterion when grading this question." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# benchmark\n", "@btime kinship($G)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Hint**: I got `@btime` output\n", "```\n", "82.144 ms (3 allocations: 7.64 MiB)\n", "```" ] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "Julia 1.7.1", "language": "julia", "name": "julia-1.7" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.7.1" }, "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 }