From d69d259dbd9d90cb4bd911a9b7273f7ab2b44d95 Mon Sep 17 00:00:00 2001 From: elvis Date: Mon, 20 Nov 2023 15:28:51 +0100 Subject: [PATCH] lesson 15/11 --- 11-15/lesson.ipynb | 259 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 259 insertions(+) create mode 100644 11-15/lesson.ipynb diff --git a/11-15/lesson.ipynb b/11-15/lesson.ipynb new file mode 100644 index 0000000..e829651 --- /dev/null +++ b/11-15/lesson.ipynb @@ -0,0 +1,259 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "0c80a449-2b98-4899-91a2-8ba36a3f73bb", + "metadata": {}, + "outputs": [], + "source": [ + "using LinearAlgebra" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "f0137444-3c22-400a-862a-8eedfbb39153", + "metadata": {}, + "outputs": [], + "source": [ + "A = [Float64(1) 1 2; 1 2 3; 3 1 4; 1 2 3+1e-8]\n", + "# very close to a singular matrix\n", + "y = A *[3;4;5] # so the solution to the lsq problem should be [3 4 5]\n", + "\n", + "F = qr(A)\n", + "U, S, V = svd(A);" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "b8736995-24a9-4488-b819-0f1afc34671b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x_qr = vcat(F.R, zeros(1, 3)) \\ ((F.Q)' * y) = [3.0000000165428684, 4.00000001654287, 4.999999983457129]\n", + "x_svd = V * (inv(diagm(S)) * (U' * y)) = [2.999999759163479, 3.999999759163469, 5.000000240836527]\n", + "x_svd2 = ((V * inv(diagm(S))) * U') * y = [2.999999669063651, 3.9999988942032676, 5.000000330936349]\n", + "x_ne = (A' * A) \\ (A' * y) = [9.864681654174836, 10.864681681143239, -1.864681661529859]\n" + ] + } + ], + "source": [ + "# some far, some close, none perfect\n", + "\n", + "@show x_qr = vcat(F.R, zeros(1,3)) \\ (F.Q'*y);\n", + "@show x_svd = V * (inv(diagm(S)) * (U' * y));\n", + "@show x_svd2 = ((V * inv(diagm(S))) * U') * y;\n", + "@show x_ne = (A'*A) \\ (A'*y);" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "fea851db-a08d-42bf-a1e1-899c63ea4b82", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "norm(A * x_qr - y) = 1.1234667099445444e-14\n", + "norm(A * x_svd - y) = 1.464821375527116e-14\n", + "norm(A * x_svd2 - y) = 2.450323676001366e-6\n", + "norm(A * x_ne - y) = 5.024432489930541e-8\n" + ] + } + ], + "source": [ + "# \"small\" errors\n", + "\n", + "@show norm(A*x_qr - y);\n", + "@show norm(A*x_svd - y);\n", + "@show norm(A*x_svd2 - y);\n", + "@show norm(A*x_ne - y);" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "06c88764-0d8c-458f-b1e3-3a5332b0df3c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2×2 Matrix{Float64}:\n", + " 2.0 1.0\n", + " 1.0 1.0" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "2×2 Matrix{Float64}:\n", + " 1.0 -1.0\n", + " -1.0 2.0" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "6.854101966249688" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "(A = [2. 1; 1 1]) |> display\n", + "display(inv(A))\n", + "cond(A) # very good" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "2ec2cc45-499b-4f36-8e30-46d834575597", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "8.55514682418072e-8" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "4.8114143668831084e-8" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "y = [2., 5]\n", + "\n", + "ytilde = y + 1e-6*randn(2)\n", + "\n", + "x = A \\ y\n", + "xtilde = A \\ ytilde\n", + "\n", + "norm(ytilde - y)/norm(y) |> display\n", + "norm(xtilde - x)/norm(x) |> display" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "29e0e48c-ec39-43a9-ae9c-16edb1ddcb97", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "400002.00000320596" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "2.767722959032909e-7" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "0.01724009150665141" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "A = [1. 1; 1 1+1e-5]\n", + "cond(A) |> display\n", + "y = [1, 1]\n", + "ytilde = y + 1e-6*randn(2)\n", + "\n", + "x = A \\ y\n", + "xtilde = A \\ ytilde\n", + "\n", + "norm(ytilde - y)/norm(y) |> display\n", + "norm(xtilde - x)/norm(x) |> display" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "eb044b9a-b0e4-4b2e-ad38-3604f5d57592", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "5.1170559300784335" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "5.693055161017966" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "A = rand(10,10)\n", + "F = svd(A)\n", + "maximum(F.S) |> display\n", + "norm(A, 2) |> display\n", + "\n", + "# they should be the same number :(" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.9.4", + "language": "julia", + "name": "julia-1.9" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.9.4" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}