{ "cells": [ { "cell_type": "markdown", "id": "53fdd565", "metadata": {}, "source": [ "# NYU Paris - Numerical Analysis" ] }, { "cell_type": "markdown", "id": "a6d022bb", "metadata": {}, "source": [ "## Final Exam\n", "\n", "- Submit this notebook on Brightspace before 6:00 PM\n", "\n", "- The exam comprises six independent exercises. In each exercise, the\n", " cells may depend on the previous cells.\n", "\n", "- To facilitate correction of your code,\n", " do not change the signatures of the functions to implement.\n", "\n", "- Run the cell below to import necessary packages and\n", " load a macro used for unit tests thoughout the notebook." ] }, { "cell_type": "code", "execution_count": null, "id": "e738e4bb", "metadata": {}, "outputs": [], "source": [ "using ForwardDiff\n", "using LinearAlgebra\n", "using Polynomials\n", "using SpecialFunctions\n", "using TestImages\n", "using LaTeXStrings\n", "using Plots\n", "\n", "macro mark(bool_expr)\n", " return :( print($bool_expr ? \"✔️\" : \"❌\"))\n", "end" ] }, { "cell_type": "markdown", "id": "0e9ebb0d", "metadata": {}, "source": [ "### 1. Floating point arithmetic\n", "\n", "Throughout this exercise,\n", "$(\\bullet)_\\beta$denotes base$\\beta$representation.\n", "True or false? (+1/0/-1)\n", "1. It holds that\n", "$\n", " (1.0111)_2 + (0.1001)_2 = (10)_2.\n", " $\n", "1. It holds that\n", "$\n", " (1000)_5 \\times (0.003)_5 = 3.\n", " $\n", "1. It holds that\n", "$\n", " (0.\\overline{4})_5 = 4\n", " $\n", "1. In base 16, all the natural numbers from 1 to 1000 can be represented using 3 digits.\n", "\n", "1. Machine multiplication according to the IEEE754 standard is a commutative operation.\n", "\n", "1. Machine addition according to the IEEE754 standard is an associative operation.\n", "\n", "1. Only finitely many real numbers can be represented exactly in the Float64 format.\n", "\n", "1. The spacing (in absolute value) between successive Float64 numbers is strictly increasing.\n", "\n", "1. In Julia, eps(Float32) is the smallest positive number representable in the Float32 format.\n", "\n", "1. In Julia, nextfloat(1f0) - 1f0 equals the machine epsilon for the Float32 format.\n", "\n", "1. Assume that$x \\in \\real$belongs to$\\mathbf F_{64}$. Then$2x \\in \\mathbf F_{64}$.\n", "\n", "1. The real number$\\sqrt{2}$can be represented exactly in the single-precision format." ] }, { "cell_type": "code", "execution_count": null, "id": "aaa7cbcc", "metadata": { "lines_to_next_cell": 2 }, "outputs": [], "source": [ "answers_q1 = zeros(Int, 12)\n", "\n", "# Use -1 for false, 1 for true\n", "answers_q1[1] = 0\n", "answers_q1[2] = 0\n", "answers_q1[3] = 0\n", "answers_q1[4] = 0\n", "answers_q1[5] = 0\n", "answers_q1[6] = 0\n", "answers_q1[7] = 0\n", "answers_q1[8] = 0\n", "answers_q1[9] = 0\n", "answers_q1[10] = 0\n", "answers_q1[11] = 0\n", "answers_q1[12] = 0" ] }, { "cell_type": "markdown", "id": "dde04a88", "metadata": {}, "source": [ "### 2. Exercise on interpolation and approximation" ] }, { "cell_type": "markdown", "id": "d7bd2802", "metadata": {}, "source": [ "1. Find the quadratic polynomial$p(x) = ax^2 + b x + c$that goes through the points$(0, 1)$,$(1, 3)$and$(2, 7),$\n", " and then plot the polynomial$p$together with the data points on the same graph." ] }, { "cell_type": "code", "execution_count": null, "id": "25bc7ff6", "metadata": {}, "outputs": [], "source": [ "# Calculate a, b and c here\n", "\n", "@mark begin p(x) = a*x^2 + b*x + c; p(0) == 1 end\n", "@mark begin p(x) = a*x^2 + b*x + c; p(1) == 3 end\n", "@mark begin p(x) = a*x^2 + b*x + c; p(2) == 7 end\n", "\n", "# Create the plot here" ] }, { "cell_type": "markdown", "id": "7b48172d", "metadata": {}, "source": [ "2. Find the function$f$of the form$f(x) = a \\sin(x) + b \\cos(x) + c \\sin(2x)$that goes through the points$(x_1, y_1)$,$(x_2, y_2)$and$(x_3, y_3)$,\n", " where the data is given below.\n", " Plot the function$f$together with the data points on the same graph." ] }, { "cell_type": "code", "execution_count": null, "id": "46a58f5a", "metadata": {}, "outputs": [], "source": [ "x = [1, 2, 3]\n", "y = [2.8313730233698577, -0.6797987415765313, -2.11828048333995]\n", "\n", "# Calculate a, b and c here\n", "\n", "@mark begin f(z) = a*sin(z) + b*cos(z) + c*sin(2z); f(x[1]) == y[1] end\n", "@mark begin f(z) = a*sin(z) + b*cos(z) + c*sin(2z); f(x[2]) == y[2] end\n", "@mark begin f(z) = a*sin(z) + b*cos(z) + c*sin(2z); f(x[3]) == y[3] end\n", "\n", "# Create the plot here" ] }, { "cell_type": "markdown", "id": "905931ce", "metadata": {}, "source": [ "3. Find the function$g$of the form$g(x) = a \\sin(x) + b \\cos(x)$\n", " such that the sum of squared residuals\n", " $$\n", " \\sum_{i=1}^{n} |g(x_i) - y_i|^2, \\qquad n = 10,\n", "$$\n", " is minimized for the data given below.\n", " Plot on the same graph the function$g$and the data points." ] }, { "cell_type": "code", "execution_count": null, "id": "bba82caa", "metadata": {}, "outputs": [], "source": [ "x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]\n", "y = [0.46618950237090034, 0.9365762460801551, 0.5907432672662498,\n", " -0.21370874419278404, -0.8869010982313607, -0.6906040605442342,\n", " 0.1784931250350807, 1.024036713535387, 0.7837248688922458,\n", " -0.1544083454499539]\n", "\n", "# Calculate a and b and create the plot here" ] }, { "cell_type": "markdown", "id": "56d0ae3b", "metadata": {}, "source": [ "### 3. Exercise on numerical integration\n", "\n", "Our goal in this exercise is to write a program to calculate integrals of the form\n", "$$\n", "I[f] := \\int_{-\\infty}^{\\infty} f(x) \\, e^{-x^2} \\, d x.\n", "$$\n", "To this end, we will use Hermite polynomials,\n", "which are orthogonal polynomials for the following inner product:\n", " $$\n", " \\langle f, g \\rangle := \\int_{-\\infty}^{\\infty} f(x) \\, g(x) \\, e^{-x^2} \\, d x\n", "$$\n", "1. Using that Hermite polynomials satisfy the recurrence relation\n", " $$\n", " H_{k + 1}(x) = 2x H_k(x) - H_k'(x), \\qquad H_0(x) = 1,\n", "$$\n", " write a function hermite(n) which returns the Hermite polynomial of degree$n$." ] }, { "cell_type": "code", "execution_count": null, "id": "4da66024", "metadata": {}, "outputs": [], "source": [ "function hermite(n)\n", " # Write your code here\n", "end\n", "\n", "@mark hermite(2) == Polynomial([-2, 0, 4])\n", "@mark hermite(3) == Polynomial([0, -12, 0, 8])\n", "@mark hermite(4) == Polynomial([12, 0, -48, 0, 16])\n", "@mark degree(hermite(10)) == 10" ] }, { "cell_type": "markdown", "id": "fde9bb30", "metadata": {}, "source": [ "2. Write a function integrate_monomial(n), which return the value of the integral\n", " $$\n", " \\int_{-\\infty}^{\\infty} x^n e^{-x^2} \\, dx =\n", " \\frac{1}{2} \\Bigl((-1)^n + 1\\Bigr) \\Gamma\\left( \\frac{n+1}{2} \\right),\n", "$$\n", " where$\\Gamma$is the Gamma function,\n", " available as gamma in the SpecialFunctions library." ] }, { "cell_type": "code", "execution_count": null, "id": "dfcc14cd", "metadata": {}, "outputs": [], "source": [ "function integrate_monomial(n)\n", " # Write your code here\n", "end\n", "\n", "@mark integrate_monomial(0) ≈ √π\n", "@mark integrate_monomial(2) ≈ √π/2\n", "@mark integrate_monomial(3) ≈ 0\n", "@mark integrate_monomial(4) ≈ 3√π/4\n", "@mark integrate_monomial(5) == 0" ] }, { "cell_type": "markdown", "id": "af088001", "metadata": {}, "source": [ "3. Write a function get_nodes_and_weights(n) which returns the nodes and weights of the Gauss-Hermite quadrature with$n$nodes.\n", "\n", " **Hints:**\n", " - Remember that the nodes of the quadrature are given by the roots of$H_{n}$.\n", " - In order to construct Lagrange polynomials, you may find it useful to use the fromroots function." ] }, { "cell_type": "code", "execution_count": null, "id": "0fc23392", "metadata": {}, "outputs": [], "source": [ "function get_nodes_and_weights(n)\n", " # Write your code here\n", " return nodes, weights\n", "end\n", "\n", "@mark length(get_nodes_and_weights(1)) == 2\n", "@mark length(get_nodes_and_weights(3)[1]) == 3\n", "@mark length(get_nodes_and_weights(5)[2]) == 5" ] }, { "cell_type": "markdown", "id": "48eb41cf", "metadata": {}, "source": [ "4. Write a function integrate_hermite(f, n), which returns an approximation of$I[f]$obtained by Gauss-Hermite integration with$n$nodes." ] }, { "cell_type": "code", "execution_count": null, "id": "dfdc0f2d", "metadata": {}, "outputs": [], "source": [ "function integrate_hermite(f, n)\n", " # Write your code here\n", "end\n", "\n", "@mark integrate_hermite(x -> 1, 10) ≈ integrate_monomial(0)\n", "@mark abs(integrate_hermite(x -> x, 10)) <= 1e-10\n", "@mark integrate_hermite(x -> x^2, 10) ≈ integrate_monomial(2)\n", "@mark integrate_hermite(x -> x^4, 10) ≈ integrate_monomial(4)\n", "@mark integrate_hermite(x -> x^6, 10) ≈ integrate_monomial(6)\n", "@mark abs(integrate_hermite(sin, 10)) <= 1e-10" ] }, { "cell_type": "markdown", "id": "308f1595", "metadata": {}, "source": [ "5. Set$n = 4$, and calculate numerically the degree of precision of the corresponding quadrature rule." ] }, { "cell_type": "code", "execution_count": null, "id": "c482007d", "metadata": {}, "outputs": [], "source": [ "# Write your code here" ] }, { "cell_type": "markdown", "id": "efd8b98a", "metadata": {}, "source": [ "6. Set$f(x) = \\cos(x)$, and plot the integration error as a function of$n$.\n", " You may use that\n", " $$\n", " \\int_{-\\infty}^{\\infty} \\cos(x) \\, e^{-x^2} \\, d x = \\frac{\\sqrt{\\pi}}{\\sqrt[4]{e}}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "a964494b", "metadata": {}, "outputs": [], "source": [ "# Write your code here" ] }, { "cell_type": "markdown", "id": "50f6b33d", "metadata": {}, "source": [ "### 4. Exercise on linear systems\n", "\n", "This exercise focuses on solving the Euler-Bernoulli beam equation in one dimension,\n", "with homogeneous Dirichlet boundary conditions:\n", "$$\n", "u''''(x) = 1, \\qquad u(0) = u'(0) = u'(1) = u(1) = 0.\n", "$$\n", "This equation models the deflection of a beam clamped at both ends,\n", "under a uniform load.\n", "A discretization of this equation on a uniform grid using the finite difference method leads to the following linear system:\n", "$$\n", "\\begin{pmatrix}\n", " 6 & -4 & 1 \\\\\n", " -4 & 6 & -4 & 1 \\\\\n", " 1 & -4 & 6 & -4 & 1 \\\\\n", " & 1 & -4 & 6 & -4 & 1 \\\\\n", " & & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots \\\\\n", " & & & 1 & -4 & 6 & -4 & 1 \\\\\n", " & & & & 1 & -4 & 6 & -4 \\\\\n", " & & & & & 1 & -4 & 6 \\\\\n", "\\end{pmatrix}\n", "\\begin{pmatrix}\n", " u_2 \\\\\n", " u_3 \\\\\n", " u_4 \\\\\n", " u_5 \\\\\n", " \\vdots \\\\\n", " u_{n-4} \\\\\n", " u_{n-3} \\\\\n", " u_{n-2}\n", "\\end{pmatrix}\n", "=\n", "h^4\n", "\\begin{pmatrix}\n", " 1 \\\\\n", " 1 \\\\\n", " 1 \\\\\n", " 1 \\\\\n", " \\vdots \\\\\n", " 1 \\\\\n", " 1 \\\\\n", " 1\n", "\\end{pmatrix},\n", "\\quad\n", "h := \\frac{1}{n},\n", "\\quad\n", "x_i := ih.\n", "$$\n", "where$h$is the spacing between the discretization points, and$(u_2, u_3, \\cdots, u_{n-3}, u_{n-2})$are the values of the unknown function$u$at the points$(x_2, x_3, \\cdots, x_{n-3}, x_{n-2})$.\n", "\n", "1. Write a function build_matrix(n), which returns the matrix of the linear system.\n", "\n", " **Hint:** the function diagm is useful here." ] }, { "cell_type": "code", "execution_count": null, "id": "0e66494d", "metadata": {}, "outputs": [], "source": [ "function build_matrix(n)\n", " # Write your code here\n", "end\n", "\n", "@mark size(build_matrix(20)) == (17, 17)\n", "@mark build_matrix(20)[5, 3] ≈ 1\n", "@mark build_matrix(20)[5, 4] ≈ -4\n", "@mark build_matrix(20)[5, 5] ≈ 6\n", "@mark build_matrix(20)[5, 6] ≈ -4\n", "@mark build_matrix(20)[5, 7] ≈ 1\n", "@mark build_matrix(20)[5, 8] ≈ 0" ] }, { "cell_type": "markdown", "id": "c67b4993", "metadata": {}, "source": [ "2. Write a function build_rhs(n), which returns the right-hand side of the linear system." ] }, { "cell_type": "code", "execution_count": null, "id": "2bf19861", "metadata": {}, "outputs": [], "source": [ "function build_rhs(n)\n", " # Write your code here\n", "end\n", "\n", "@mark length(build_rhs(20)) == 17\n", "@mark build_rhs(20)[end] == 1 / 20^4\n", "@mark build_rhs(20)[1] == 1 / 20^4" ] }, { "cell_type": "markdown", "id": "e5de6676", "metadata": {}, "source": [ "3. Write a function forward_substitution(L, y)\n", " that solves the linear system$\\mathsf L \\mathbf x = \\mathbf y$,\n", " in the case where L is a lower-triangular matrix,\n", " by using forward substitution." ] }, { "cell_type": "code", "execution_count": null, "id": "a789548e", "metadata": {}, "outputs": [], "source": [ "function forward_substitution(L, y)\n", " # Write your code here\n", "end\n", "\n", "@mark begin n = 10; L = [1 0; 2 1]; forward_substitution(L, [1; 0]) ≈ [1; -2] end\n", "@mark begin n = 10; L = LowerTriangular(ones(n, n)); y = fill(2, n); forward_substitution(L, L*y) ≈ y end\n", "@mark begin n = 10; L = LowerTriangular(randn(n, n)); y = randn(n); forward_substitution(L, L*y) ≈ y end\n", "@mark begin n = 20; L = LowerTriangular(2ones(n, n)); y = rand(n); forward_substitution(L, L*y) ≈ y end" ] }, { "cell_type": "markdown", "id": "10eeaf06", "metadata": {}, "source": [ "4. The successive over-relaxation method is a splitting method for solving linear systems of the form$\\mathsf A \\mathbf x = \\mathbf b$.\n", " It is based on the iteration\n", " $$\n", " \\mathsf M \\mathbf x_{k+1} = \\mathsf N \\mathbf x_{k} + \\mathbf{b}, \\qquad \\text{where} \\quad \\mathsf M = \\frac{\\mathsf D}{\\omega} + \\mathsf L \\quad \\text{and} \\quad \\mathsf N = \\mathsf M - \\mathsf A.\n", " \\tag{SOR}\n", "$$\n", " \n", " Here$\\omega \\in (0, 2)$is a parameter,\n", "$\\mathsf D$is diagonal matrix containing the diagonal of$\\mathsf A$,\n", " and$\\mathsf L$is a strictly lower triangular matrix containing the strict lower triangular part of$\\mathsf A$,\n", " not including the diagonal.\n", " Write a function sor(A, b, ω) that\n", " implements this iteration, without using Julia's \\ and inv functions.\n", " Initialize the iteration with$\\mathbf x_0 = \\mathbf 0$,\n", " and use as stopping criterion that\n", " $$\n", " \\lVert \\mathsf A \\mathbf x - \\mathbf b \\rVert\n", " \\leq \\varepsilon \\lVert \\mathbf b \\rVert,\n", " \\qquad \\varepsilon := 10^{-10}.\n", "$$\n", " If after maxiter iterations,\n", " a solution that satisfies this stopping criterion has not been found,\n", " return nothing.\n", "\n", " **Hints**:\n", " - The functions Diagonal and LowerTriangular are useful to construct the matrices$\\mathsf D$and$\\mathsf L$.\n", " - In order to solve (SOR) at each iteration,\n", " use the function forward_substitution you wrote above." ] }, { "cell_type": "code", "execution_count": null, "id": "735e42a1", "metadata": {}, "outputs": [], "source": [ "function sor(A, b; ω = 1, ε = 1e-10, maxiter = 10^5)\n", " # Write your code here\n", "end\n", "\n", "# Test code\n", "n = 20\n", "X = qr(randn(n, n)).Q\n", "A = X*Diagonal(rand(n) .+ 5)*X'\n", "b = randn(n)\n", "@mark begin ε = 1e-10; sol = sor(A, b; ω = 1.5, ε = ε); norm(A*sol - b) ≤ ε*norm(b) end\n", "@mark begin ε = 1e-10; sol = sor(A, b; ω = 1.5, ε = ε); norm(A*sol - b) ≥ 1e-15*norm(b) end\n", "@mark begin ε = 1e-12; sol = sor(A, b; ω = 1.5, ε = ε); norm(A*sol - b) ≤ ε*norm(b) end\n", "@mark begin ε = 1e-12; sor(A, b; ω = 2, ε = ε) == nothing end\n", "@mark begin ε = 1e-12; sor(A, b; ω = 1, ε = ε) ≈ A\\b end" ] }, { "cell_type": "markdown", "id": "43e63cfc", "metadata": {}, "source": [ "5. Use the relaxation method implemented in the previous item\n", " with parameters$\\omega = 1.5$and$\\varepsilon = 10^{-8}$to solve the linear system of this exercise with$n = 40$.\n", " Compare on a graph the solution you obatined with the exact solution,\n", " which is given by\n", " $$u(x) = \\frac{1}{24} x^2(x - 1)^2.$$" ] }, { "cell_type": "code", "execution_count": null, "id": "1f5ba663", "metadata": { "lines_to_next_cell": 2 }, "outputs": [], "source": [ "u(x) = 1/24 * x^2 * (x - 1)^2\n", "# Write your code here" ] }, { "cell_type": "markdown", "id": "badd118f", "metadata": {}, "source": [ "### 5. Exercise on nonlinear equations\n", "\n", "We wish to find all the solutions to the following transcendental equation for$x \\in [-5, -5]$.\n", "$$\n", "x = 5 \\sin(\\pi x)\n", "\\tag{1}\n", "$$\n", "\n", "1. Plot the functions$x \\mapsto x$and$x \\mapsto 5 \\sin(\\pi x)$on the same graph,\n", " for$x$in the range$[-5, 5]$,\n", " and count the number of solutions" ] }, { "cell_type": "code", "execution_count": null, "id": "a8060fe7", "metadata": {}, "outputs": [], "source": [ "# Write your code here" ] }, { "cell_type": "markdown", "id": "7803e926", "metadata": {}, "source": [ "2. Using the bisection method, calculate precisely the only solution in the interval$[\\frac{1}{2}, 1]$." ] }, { "cell_type": "code", "execution_count": null, "id": "969f0622", "metadata": {}, "outputs": [], "source": [ "function bisection(f, a, b, ε = 1e-10)\n", " # Write your code here\n", "end\n", "\n", "# Calculate solution here" ] }, { "cell_type": "markdown", "id": "6e92a82a", "metadata": {}, "source": [ "3. Write a function newton_raphson(f, x₀; maxiter = 100, ε = 1e-12) that returns the result of the Newton-Raphson iteration applied to the equation$f(x) = 0$and initialized at x₀,\n", " or nothing if a solution is not found after maxiter iterations.\n", " Use the ForwardDiff library to compute derivatives,\n", " and use the following stopping criterion\n", " $$\n", " |f(x_k)| ≤ \\varepsilon.\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "a2f62f9d", "metadata": {}, "outputs": [], "source": [ "function newton_raphson(f, x₀; maxiter = 100, ε = 1e-12)\n", " # Write your code here\n", "end\n", "\n", "@mark newton_raphson(x -> x^2 - 2, 1) ≈ √2\n", "@mark newton_raphson(x -> x^2 - 2, -1) ≈ -√2\n", "@mark newton_raphson(x -> x^3 - 2, 1) ≈ cbrt(2)\n", "@mark newton_raphson(x -> x^2 + 1, 1.5) == nothing\n", "@mark newton_raphson(x -> x^2 + 1, 1) == nothing" ] }, { "cell_type": "markdown", "id": "2219fc5f", "metadata": {}, "source": [ "4. Using the Newton-Raphson method you implemented,\n", " calculate all the solutions to the transcendental equation (1)." ] }, { "cell_type": "code", "execution_count": null, "id": "ce44d968", "metadata": {}, "outputs": [], "source": [ "# Write your code here." ] }, { "cell_type": "markdown", "id": "d4bcb668", "metadata": {}, "source": [ "5. We now consider another approach,\n", " called the *secant method*.\n", " An iteration of this method is given by\n", " $$\n", " x_{k+2} = \\frac{x_{k} f(x_{k+1}) - x_{k+1}f(x_k)}{f(x_{k+1}) - f(x_{k})}.\n", "$$\n", " Write a function secant(f, x₀, x₁; maxiter = 100, ε = 1e-12)\n", " that returns the result of the secant iteration applied to the equation$f(x) = 0$, and initialized with x₀ and x₁,\n", " or nothing if a solution is not found after maxiter iterations.\n", " Use the same stopping criterion as above." ] }, { "cell_type": "code", "execution_count": null, "id": "8819b0c1", "metadata": {}, "outputs": [], "source": [ "function secant(f, x₀, x₁; maxiter = 100, ε = 1e-12)\n", " # Write your code here\n", "end\n", "\n", "@mark secant(x -> x^2 - 2, 1., 2.) ≈ √2\n", "@mark secant(x -> x^2 - 2, -1., -2.) ≈ -√2\n", "@mark secant(x -> x^2 + 1, 1., 2.) == nothing" ] }, { "cell_type": "markdown", "id": "48435545", "metadata": {}, "source": [ "6. Use the secant method you implemented to solve the transcendental equation\n", "$$x = e^{-x}$$" ] }, { "cell_type": "code", "execution_count": null, "id": "bb25e735", "metadata": {}, "outputs": [], "source": [ "# Write your code here." ] }, { "cell_type": "markdown", "id": "f6721f3b", "metadata": {}, "source": [ "### 6. Exercise on eigenvalue problems" ] }, { "cell_type": "markdown", "id": "cb8d2606", "metadata": {}, "source": [ "Our goal in this exercise is to implement,\n", "without using the functions eigen and svd from the LinearAlgebra library,\n", "a simple algorithm for calculating the dominant *singular values* and associated *singular vectors* of a matrix,\n", "and to use this algorithm for image compression.\n", "Recall that any matrix$\\mathsf A \\in \\mathbf R^{n \\times n}$admits a *singular value decomposition* of the form\n", "$$\n", "\\mathsf A = \\mathsf U \\mathsf \\Sigma \\mathsf V^\\top,\n", "$$\n", "where$\\mathsf U \\in \\mathbf R^{n \\times n}$and$\\mathsf V \\in \\mathbf R^{n \\times n}$are orthogonol matrices,\n", "and$\\mathsf \\Sigma \\in \\mathbf R^{n \\times n}$is a diagonal matrix.\n", "The diagonal entries of$\\mathsf \\Sigma$are called *singular value*,\n", "while the columns of$\\mathsf U$and$\\mathsf V$are the associated left and right *singular vectors*, respectively.\n", "We restrict our attention to square matrices, for simplicity.\n", "\n", "1. Write a function myEigen(B, p, niter) which returns the p dominant eigenvalues,\n", " together with the associated eigenvectors,\n", " of a symmetric matrix B.\n", " To this end, implement the simultaneous iteration seen in class.\n", " You can use the myQR function given below in order to calculate reduced QR decompositions." ] }, { "cell_type": "code", "execution_count": null, "id": "c1dce274", "metadata": {}, "outputs": [], "source": [ "function myQR(B)\n", " Q, R = qr(B)\n", " return Matrix(Q), R\n", "end\n", "\n", "function myEigen(B, p, niter)\n", " # Write your code here\n", "end\n", "\n", "@mark begin n = 2; A = randn(n, n); myEigen(A'A, n, 100)[1] ≈ reverse(eigen(A'A).values) end\n", "@mark begin n = 4; A = randn(n, n); myEigen(A'A, n, 100)[1] ≈ reverse(eigen(A'A).values) end\n", "@mark begin A = randn(5, 5); q, r = qr(A); B = q*Diagonal(1:5)*q'; myEigen(B, 3, 100)[1] ≈ [5; 4; 3] end" ] }, { "cell_type": "markdown", "id": "6ec0399d", "metadata": {}, "source": [ "2. Write a function mySVD(B, p, niter)\n", " that returns the p dominant singular values of square matrix B (in a vector σs),\n", " together with the associated left and right singular vectors (in matrices Up and Vp).\n", " To this end,\n", " notice that\n", " $$\n", " \\mathsf A \\mathsf A^\\top = \\mathsf U \\mathsf \\Sigma^2 \\mathsf U^\\top, \\qquad\n", " \\mathsf A^\\top \\mathsf A = \\mathsf V \\mathsf \\Sigma^2 \\mathsf V^\\top.\n", "$$\n", " Therefore, the left singular vectors of$\\mathsf A$\n", " are the eigenvectors of$\\mathsf A \\mathsf A^\\top$,\n", " while the right singular vectors of$\\mathsf A$\n", " are the eigenvectors of$\\mathsf A^\\top \\mathsf A$,\n", "\n", " **Hints**:\n", "\n", " - In order to calculate singular vectors,\n", " apply the function myEigen you wrote above to the matrices$\\mathsf A \\mathsf A^\\top$\n", " (in order to obtain left singular vectors),\n", " and$\\mathsf A^\\top \\mathsf A$(in order to obtain right singular vectors).\n", "\n", " - Once you have calculated the left and right singular vectors\n", " associated with the p dominant singular values,\n", " the singular values themselves can be obtained by extracting the diagonal from the matrix\n", " $$\n", " \\Sigma_p = \\mathsf U_p^\\top \\mathsf B \\mathsf V_p.\n", "$$\n", " Here$\\mathsf U_p$and$\\mathsf V_p$are the matrices containing as columns the left and right singular vectors associated with the p dominant singular values,\n", " respectively." ] }, { "cell_type": "code", "execution_count": null, "id": "991a28c2", "metadata": {}, "outputs": [], "source": [ "function mySVD(B, p, niter)\n", " # Write your code here\n", "end" ] }, { "cell_type": "markdown", "id": "1dfe8108", "metadata": {}, "source": [ "3. Write a function to test your implementation of mySVD.\n", " Proceed as follows:\n", " - Create a 10 x 10 matrix B with random entries using the randn function.\n", " - Calculate its full singular value decomposition using mySVD.\n", " - Check that U and V are orthogonal,\n", " and that U*Diagonal(σs)*V' is close to B." ] }, { "cell_type": "code", "execution_count": null, "id": "81b070a6", "metadata": {}, "outputs": [], "source": [ "# Write you code here" ] }, { "cell_type": "markdown", "id": "e0113e32", "metadata": {}, "source": [ "4. The singular value decomposition is very useful for compressing matrices.\n", " The idea of matrix compression based on SVD is the following:\n", " given$p \\leqslant n$, a matrix$B \\in \\mathbf R^{n\\times n}$\n", " can be approximated by\n", " $$\n", " \\widetilde {\\mathsf B} := \\mathsf U_p \\mathsf \\Sigma_p \\mathsf V_p^\\top,\n", "$$\n", " where$\\Sigma_p \\in \\mathbf R^{p \\times p}$is a diagonal matrix containing the$p$dominant singular values of$\\mathsf B$on its diagonal,\n", " and where$\\mathsf U_p \\in \\mathbf R^{n \\times p}$and$\\mathsf V_p \\in \\mathbf R^{n \\times p}$are rectangular matrices containing the associated left and right singular vectors,\n", " respectively.\n", "\n", " Since a grayscale image can be represented by a matrix containing the intensity values of all the pixels,\n", " this approach for compressing matrices can be used for compressing grayscale images.\n", " Use this method, i.e. calculate$\\widetilde {\\mathsf B}$, for$p \\in \\{5, 10, 20, 30\\}$,\n", " in order to compress the image fabio_gray_256 given below,\n", " and plot the compressed image for these values of p." ] }, { "cell_type": "markdown", "id": "01a449fc", "metadata": {}, "source": [ " **Remarks**:\n", " - (For information only) In practice, instead of storing the full matrix$\\widetilde {\\mathsf B}$, which contains$n^2$entries,\n", " we can store only the matrices$\\mathsf U_p$,$\\mathsf \\Sigma_p$and$\\mathsf V_p$,\n", " which together contain only$2np + p^2$entries.\n", " If$p \\ll n$,\n", " then the memory required to store these matrices is much smaller than the memory required to store the initial matrix$\\mathsf B\$.\n", "\n", " - A function for drawing images based on the matrix of pixel intensity values is provided below." ] }, { "cell_type": "code", "execution_count": null, "id": "612d3d33", "metadata": {}, "outputs": [], "source": [ "# Do not change the code in this cell\n", "A = testimage(\"fabio_gray_256\")\n", "\n", "# Convert image to matrix of Float64\n", "M = Float64.(A)\n", "\n", "# Function to plot a grayscale image from the matrix\n", "# containing the intensity values of all the pixels\n", "function plot_matrix(B)\n", " plot(Gray.(B), ticks=false, showaxis=false)\n", "end\n", "\n", "plot_matrix(M)" ] }, { "cell_type": "code", "execution_count": null, "id": "56bda9b6", "metadata": {}, "outputs": [], "source": [ "p = 5\n", "# Write you code here" ] }, { "cell_type": "code", "execution_count": null, "id": "5801459e", "metadata": {}, "outputs": [], "source": [ "p = 10\n", "# Write you code here" ] }, { "cell_type": "code", "execution_count": null, "id": "b028396c", "metadata": {}, "outputs": [], "source": [ "p = 20\n", "# Write you code here" ] }, { "cell_type": "code", "execution_count": null, "id": "a3d591f0", "metadata": {}, "outputs": [], "source": [ "p = 30\n", "# Write you code here" ] } ], "metadata": { "jupytext": { "encoding": "# -*- coding: utf-8 -*-", "main_language": "julia" } }, "nbformat": 4, "nbformat_minor": 5 }