\n",
"```\n",
"\n",
"```{raw} latex\n",
"%%start demo%%\n",
"```\n",
"\n",
"```{index} ! Julia; fill, Julia; diagm, ! Julia; diag\n",
"```\n",
"\n",
"::::{panels}\n",
":column: col-7 left-side\n",
":card: border-0 shadow-none\n",
"```{raw} latex\n",
"\\begin{minipage}[t]{0.5\\textwidth}\n",
"```\n",
"Here is a small tridiagonal matrix. Note that there is one fewer element on the super- and subdiagonals than on the main diagonal.\n",
"\n",
"```{raw} latex\n",
"\\end{minipage}\\hfill\n",
"```\n",
"---\n",
":column: col-5 right-side\n",
":card: shadow-none comment\n",
"```{raw} latex\n",
"\\begin{minipage}[t]{0.4\\textwidth}\\begin{mdframed}[default]\\small\n",
"```\n",
"Use `fill` to create an array of a given size, with each element equal to a provided value.\n",
"```{raw} latex\n",
"\\end{mdframed}\\end{minipage}\n",
"```\n",
"::::"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "27f04d79",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6×6 Matrix{Int64}:\n",
" 2 -1 0 0 0 0\n",
" 4 2 -1 0 0 0\n",
" 0 3 0 -1 0 0\n",
" 0 0 2 2 -1 0\n",
" 0 0 0 1 1 -1\n",
" 0 0 0 0 0 2"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = diagm( -1=>[4,3,2,1,0], 0=>[2,2,0,2,1,2], 1=>fill(-1,5) )"
]
},
{
"cell_type": "markdown",
"id": "57cc5c59",
"metadata": {},
"source": [
"```{index} ! Julia; diag\n",
"```\n",
"\n",
"::::{panels}\n",
":column: col-7 left-side\n",
":card: border-0 shadow-none\n",
"```{raw} latex\n",
"\\begin{minipage}[t]{0.5\\textwidth}\n",
"```\n",
"We can extract the elements on any diagonal using the `diag` function. The main or central diagonal is numbered zero, above and to the right of that is positive, and below and to the left is negative.\n",
"```{raw} latex\n",
"\\end{minipage}\\hfill\n",
"```\n",
"---\n",
":column: col-5 right-side\n",
":card: shadow-none comment\n",
"```{raw} latex\n",
"\\begin{minipage}[t]{0.4\\textwidth}\\begin{mdframed}[default]\\small\n",
"```\n",
"The `diag` function extracts the elements from a specified diagonal of a matrix.\n",
"```{raw} latex\n",
"\\end{mdframed}\\end{minipage}\n",
"```\n",
"::::"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "7fd1086e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"diag_main = diag(A) = [2, 2, 0, 2, 1, 2]\n",
"diag_minusone = diag(A, -1) = [4, 3, 2, 1, 0]\n"
]
}
],
"source": [
"@show diag_main = diag(A);\n",
"@show diag_minusone = diag(A,-1);"
]
},
{
"cell_type": "markdown",
"id": "5057af77",
"metadata": {},
"source": [
"The lower and upper bandwidths of $\\mathbf{A}$ are repeated in the factors from the unpivoted LU factorization."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "3f0c6488",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6×6 LowerTriangular{Float64, Matrix{Float64}}:\n",
" 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ \n",
" 2.0 1.0 ⋅ ⋅ ⋅ ⋅ \n",
" 0.0 0.75 1.0 ⋅ ⋅ ⋅ \n",
" 0.0 0.0 2.66667 1.0 ⋅ ⋅ \n",
" 0.0 0.0 0.0 0.214286 1.0 ⋅ \n",
" 0.0 0.0 0.0 0.0 0.0 1.0"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"L,U = FNC.lufact(A)\n",
"L"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "91bc011d",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6×6 UpperTriangular{Float64, Matrix{Float64}}:\n",
" 2.0 -1.0 0.0 0.0 0.0 0.0\n",
" ⋅ 4.0 -1.0 0.0 0.0 0.0\n",
" ⋅ ⋅ 0.75 -1.0 0.0 0.0\n",
" ⋅ ⋅ ⋅ 4.66667 -1.0 0.0\n",
" ⋅ ⋅ ⋅ ⋅ 1.21429 -1.0\n",
" ⋅ ⋅ ⋅ ⋅ ⋅ 2.0"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"U"
]
},
{
"cell_type": "markdown",
"id": "d465720c",
"metadata": {},
"source": [
"```{raw} html\n",
"

\n",
"```\n",
"\n",
"```{raw} latex\n",
"%%end demo%%\n",
"```\n",
"\n",
"```{index} pivoting\n",
"```\n",
"\n",
"If row pivoting is not used, the $\\mathbf{L}$ and $\\mathbf{U}$ factors preserve the lower and upper bandwidths of $\\mathbf{A}$. This fact implies computational savings in both the factorization and the triangular substitutions because the zeros appear predictably and we can skip multiplication and addition with them. \n",
"\n",
"```{proof:observation}\n",
"The number of flops needed by LU factorization without pivoting is $O(b_u b_\\ell n)$ when the upper and lower bandwidths are $b_u$ and $b_\\ell$.\n",
"```\n",
"\n",
"```{index} sparse matrix\n",
"```\n",
"\n",
"In order to exploit the savings offered by sparsity, we would need to make modifications to {numref}`Function %s \n",
"```\n",
"\n",
"```{raw} latex\n",
"%%start demo%%\n",
"```\n",
"\n",
"If we use an ordinary or dense matrix, then there's no way to exploit a banded structure such as tridiagonality."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "cb5f124c",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 8.652754 seconds (5 allocations: 763.016 MiB, 0.23% gc time)\n"
]
}
],
"source": [
"n = 10000\n",
"A = diagm(0=>1:n, 1=>n-1:-1:1, -1=>ones(n-1))\n",
"lu(rand(3,3)) # force compilation\n",
"@time lu(A);"
]
},
{
"cell_type": "markdown",
"id": "63fd79be",
"metadata": {},
"source": [
"If instead we construct a proper sparse matrix, though, the speedup can be dramatic."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "66659b65",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.005833 seconds (65 allocations: 9.766 MiB)\n"
]
}
],
"source": [
"A = spdiagm(0=>1:n, 1=>n-1:-1:1, -1=>ones(n-1))\n",
"lu(A); # compile for sparse case\n",
"@time lu(A);"
]
},
{
"cell_type": "markdown",
"id": "9d24d25d",
"metadata": {},
"source": [
"You can also see above that far less memory was used in the sparse case.\n",
"```{raw} html\n",
"

\n",
"```\n",
"\n",
"```{raw} latex\n",
"%%end demo%%\n",
"```\n",
"\n",
"## Symmetric matrices\n",
"\n",
"```{index} symmetric matrix\n",
"```\n",
"\n",
"::::{proof:definition} Symmetric matrix\n",
"A square matrix $\\mathbf{A}$ satisfying $\\mathbf{A}^T = \\mathbf{A}$ is called **symmetric**.\n",
"::::\n",
"\n",
"Symmetric matrices arise frequently in applications because many types of interactions, such as gravitation and social-network befriending, are inherently symmetric. Symmetry in linear algebra simplifies many properties and algorithms. As a rule of thumb, if your matrix has symmetry, you want to exploit and preserve it. \n",
"\n",
"In $\\mathbf{A}=\\mathbf{L}\\mathbf{U}$ we arbitrarily required the diagonal elements of $\\mathbf{L}$, but not $\\mathbf{U}$, to be one. That breaks symmetry, so we need to modify the goal to\n",
"\n",
"$$\n",
"\\mathbf{A}=\\mathbf{L}\\mathbf{D}\\mathbf{L}^T,\n",
"$$\n",
"\n",
"where $\\mathbf{L}$ is unit lower triangular and $\\mathbf{D}$ is diagonal. To find an algorithm for this factorization, we begin by generalizing {eq}`matrixouter` a bit without furnishing proof.\n",
"\n",
"::::{proof:theorem} Linear combination of outer products\n",
"Let $\\mathbf{D}$ be an $n\\times n$ diagonal matrix with diagonal elements $d_1,d_2,\\ldots,d_n$, and suppose $\\mathbf{A}$ and $\\mathbf{B}$ are $n\\times n$ as well. Write the columns of $\\mathbf{A}$ as $\\mathbf{a}_1,\\dots,\\mathbf{a}_n$ and the rows of $\\mathbf{B}$ as $\\mathbf{b}_1^T,\\dots,\\mathbf{b}_n^T$. Then\n",
"\n",
"```{math}\n",
":label: matrixouter3\n",
"\\mathbf{A}\\mathbf{D}\\mathbf{B} = \\sum_{k=1}^n d_k \\mathbf{a}_k \\mathbf{b}_k^T.\n",
"```\n",
"::::\n",
"\n",
"Let's derive the LDL$^T$ factorization for a small example.\n",
"\n",
"(demo-structure-symm)=\n",
"```{proof:demo}\n",
"```\n",
"\n",
"```{raw} html\n",
"\n",
"```\n",
"\n",
"```{raw} latex\n",
"%%start demo%%\n",
"```\n",
"\n",
"We begin with a symmetric $\\mathbf{A}$."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "810daed0",
"metadata": {},
"outputs": [],
"source": [
"A₁ = [ 2 4 4 2\n",
" 4 5 8 -5\n",
" 4 8 6 2\n",
" 2 -5 2 -26 ];"
]
},
{
"cell_type": "markdown",
"id": "fe8c2397",
"metadata": {},
"source": [
"We won't use pivoting, so the pivot element is at position (1,1). This will become the first element on the diagonal of $\\mathbf{D}$. Then we divide by that pivot to get the first column of $\\mathbf{L}$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "fc40ac61",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Matrix{Float64}:\n",
" 0.0 0.0 0.0 0.0\n",
" 0.0 -3.0 0.0 -9.0\n",
" 0.0 0.0 -2.0 -2.0\n",
" 0.0 -9.0 -2.0 -28.0"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"L = diagm(ones(4))\n",
"d = zeros(4)\n",
"d[1] = A₁[1,1]\n",
"L[:,1] = A₁[:,1]/d[1]\n",
"A₂ = A₁ - d[1]*L[:,1]*L[:,1]'"
]
},
{
"cell_type": "markdown",
"id": "e1d8b805",
"metadata": {},
"source": [
"We are now set up the same way for the submatrix in rows and columns 2–4."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "06760e83",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Matrix{Float64}:\n",
" 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 -2.0 -2.0\n",
" 0.0 0.0 -2.0 -1.0"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"d[2] = A₂[2,2]\n",
"L[:,2] = A₂[:,2]/d[2]\n",
"A₃ = A₂ - d[2]*L[:,2]*L[:,2]'"
]
},
{
"cell_type": "markdown",
"id": "3ec93d74",
"metadata": {},
"source": [
"We continue working our way down the diagonal."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "a79cb741",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"d = [2.0, -3.0, -2.0, 1.0]\n"
]
},
{
"data": {
"text/plain": [
"4×4 Matrix{Float64}:\n",
" 1.0 -0.0 -0.0 0.0\n",
" 2.0 1.0 -0.0 0.0\n",
" 2.0 -0.0 1.0 0.0\n",
" 1.0 3.0 1.0 1.0"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"d[3] = A₃[3,3]\n",
"L[:,3] = A₃[:,3]/d[3]\n",
"A₄ = A₃ - d[3]*L[:,3]*L[:,3]'\n",
"d[4] = A₄[4,4]\n",
"@show d;\n",
"L"
]
},
{
"cell_type": "markdown",
"id": "c65f2720",
"metadata": {},
"source": [
"We have arrived at the targeted factorization."
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "a40f3e48",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Matrix{Float64}:\n",
" 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A₁ - L*diagm(d)*L'"
]
},
{
"cell_type": "markdown",
"id": "371acc79",
"metadata": {},
"source": [
"```{raw} html\n",
"

\n",
"```\n",
"\n",
"```{raw} latex\n",
"%%end demo%%\n",
"```\n",
"\n",
"In practice we don't actually have to carry out any arithmetic in the upper triangle of $\\mathbf{A}$ as we work, since the operations are always the mirror image of those in the lower triangle. As a result, it can be shown that LDL$^T$ factorization takes about half as much work as the standard LU.\n",
"\n",
"::::{proof:observation}\n",
"LDL$^T$ factorization on an $n \\times n$ symmetric matrix, when successful, takes $\\sim \\frac{1}{3}n^3$ flops.\n",
"::::\n",
"\n",
"Just as pivoting is necessary to stabilize LU factorization, the LDL$^T$ factorization without pivoting may be unstable or even fail to exist. We won't go into the details, because our interest is in specializing the factorization to matrices that also possess another important property.\n",
"\n",
"(sec-SPD)=\n",
"\n",
"## Symmetric positive definite matrices\n",
"\n",
"Suppose that $\\mathbf{A}$ is $n\\times n$ and $\\mathbf{x}\\in\\mathbb{R}^n$. Observe that $\\mathbf{x}^T\\mathbf{A}\\mathbf{x}$ is the product of $1\\times n$, $n\\times n$, and $n\\times 1$ matrices, so it is a scalar, sometimes referred to as a **quadratic form**. It can be expressed as\n",
"\n",
"```{math}\n",
" :label: quadratic-form\n",
" \\mathbf{x}^T\\mathbf{A}\\mathbf{x} = \\sum_{i=1}^n \\sum_{j=1}^n A_{ij}x_ix_j.\n",
"```\n",
"\n",
"```{index} ! symmetric positive definite matrix\n",
"```\n",
"\n",
"```{index} see: SPD matrix; symmetric positive definite matrix\n",
"```\n",
"\n",
"::::{proof:definition} Symmetric positive definite matrix\n",
"A real $n\\times n$ matrix $\\mathbf{A}$ is called a **symmetric positive definite matrix** (or SPD matrix) if it is symmetric and, for all nonzero $\\mathbf{x}\\in\\mathbb{R}^n$,\n",
"\n",
"```{math}\n",
" :label: SPD-def\n",
" \\mathbf{x}^T \\mathbf{A} \\mathbf{x} > 0.\n",
"```\n",
"::::\n",
"\n",
"The definiteness property is usually difficult to check directly from the definition. There are some equivalent conditions, though. For instance, a symmetric matrix is positive definite if and only if its eigenvalues are all real positive numbers. SPD matrices have important properties and appear in applications in which the definiteness is known for theoretical reasons.\n",
"\n",
"Let us consider what definiteness means to the LDL$^T$ factorization. We compute\n",
"\n",
"```{math}\n",
" 0 < \\mathbf{x}^T\\mathbf{A}\\mathbf{x} = \\mathbf{x}^T \\mathbf{L} \\mathbf{D} \\mathbf{L}^T \\mathbf{x} = \\mathbf{z}^T \\mathbf{D} \\mathbf{z},\n",
"```\n",
"\n",
"where $\\mathbf{z}=\\mathbf{L}^T \\mathbf{x}$. Note that since $\\mathbf{L}$ is unit lower triangular, it is nonsingular, so $\\mathbf{x}=\\mathbf{L}^{-T}\\mathbf{z}$. By taking $\\mathbf{z}=\\mathbf{e}_k$ for $k=1,\\ldots,n$, we can read the equalities from right to left and conclude that $D_{kk}>0$ for all $k$. That permits us to write a kind of square root formula:[^sqrt]\n",
"\n",
"[^sqrt]: Except for this diagonal, positive definite case, it's not trivial to define the square root of a matrix, so don't generalize the notation used here.\n",
"\n",
"```{math}\n",
" :label: diag-sqrt\n",
" \\mathbf{D} =\n",
" \\begin{bmatrix}\n",
" D_{11} & & & \\\\\n",
" & D_{22} & & \\\\\n",
" & & \\ddots & \\\\\n",
" & & & D_{nn}\n",
" \\end{bmatrix}\n",
"= \\begin{bmatrix}\n",
" \\sqrt{D_{11}} & & & \\\\\n",
" & \\sqrt{D_{22}} & & \\\\\n",
" & & \\ddots & \\\\\n",
" & & & \\sqrt{D_{nn}}\n",
" \\end{bmatrix}^{\\,2}\n",
"= \\bigl( \\mathbf{D}^{1/2} \\bigr)^2.\n",
"```\n",
"\n",
"Now we have $\\mathbf{A}=\\mathbf{L}\\mathbf{D}^{1/2}\\mathbf{D}^{1/2}\\mathbf{L}^T= \\mathbf{R}^T \\mathbf{R}$, where $\\mathbf{R} =\\mathbf{D}^{1/2}\\mathbf{L}^T$ is an upper triangular matrix whose diagonal entries are positive.\n",
"\n",
"```{index} ! matrix factorization; Cholesky\n",
"```\n",
"\n",
"```{index} pivoting\n",
"```\n",
"\n",
"::::{proof:theorem} Cholesky factorization\n",
"Any SPD matrix $\\mathbf{A}$ may be factored as \n",
"\n",
"$$\n",
"\\mathbf{A} = \\mathbf{R}^T \\mathbf{R},\n",
"$$\n",
"\n",
"where $\\mathbf{R}$ is an upper triangular matrix with positive diagonal elements. This is called the **Cholesky factorization**.\n",
"::::\n",
"\n",
"While the unpivoted LDL$^T$ factorization is not stable and not even always possible, in the SPD case one can prove that pivoting is not necessary for the existence nor the stability of the Cholesky factorization. \n",
"\n",
"::::{proof:observation}\n",
"Cholesky factorization of an $n \\times n$ SPD matrix takes $\\sim \\frac{1}{3}n^3$ flops.\n",
"::::\n",
"\n",
"The speed and stability of the Cholesky factorization make it the top choice for solving linear systems with SPD matrices. As a side benefit, the Cholesky algorithm fails (in the form of an imaginary square root or division by zero) if and only if the matrix $\\mathbf{A}$ is not positive definite. This is often the best way to test the definiteness of a symmetric matrix about which nothing else is known.\n",
"\n",
"(demo-structure-cholesky)=\n",
"```{proof:demo}\n",
"```\n",
"\n",
"```{raw} html\n",
"\n",
"```\n",
"\n",
"```{raw} latex\n",
"%%start demo%%\n",
"```\n",
"\n",
"\n",
"A randomly chosen matrix is extremely unlikely to be symmetric. However, there is a simple way to symmetrize one."
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "ca0e68f2",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Matrix{Float64}:\n",
" 8.0 14.0 9.0 8.0\n",
" 14.0 12.0 8.0 10.0\n",
" 9.0 8.0 6.0 14.0\n",
" 8.0 10.0 14.0 14.0"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = rand(1.0:9.0,4,4)\n",
"B = A + A'"
]
},
{
"cell_type": "markdown",
"id": "7b436e4e",
"metadata": {},
"source": [
"```{index} ! Julia; cholesky\n",
"```\n",
"\n",
"::::{panels}\n",
":column: col-7 left-side\n",
":card: border-0 shadow-none\n",
"```{raw} latex\n",
"\\begin{minipage}[t]{0.5\\textwidth}\n",
"```\n",
"Similarly, a random symmetric matrix is unlikely to be positive definite. The Cholesky algorithm always detects a non-PD matrix by quitting with an error.\n",
"\n",
"```{raw} latex\n",
"\\end{minipage}\\hfill\n",
"```\n",
"---\n",
":column: col-5 right-side\n",
":card: shadow-none comment\n",
"```{raw} latex\n",
"\\begin{minipage}[t]{0.4\\textwidth}\\begin{mdframed}[default]\\small\n",
"```\n",
"The `cholesky` function computes a Cholesky factorization if possible, or throws an error for a non-positive-definite matrix. However, it does *not* check for symmetry.\n",
"```{raw} latex\n",
"\\end{mdframed}\\end{minipage}\n",
"```\n",
"::::"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "4c720103",
"metadata": {},
"outputs": [
{
"ename": "LoadError",
"evalue": "PosDefException: matrix is not positive definite; Cholesky factorization failed.",
"output_type": "error",
"traceback": [
"PosDefException: matrix is not positive definite; Cholesky factorization failed.",
"",
"Stacktrace:",
" [1] checkpositivedefinite",
" @ /Applications/julia-rosetta/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/factorization.jl:18 [inlined]",
" [2] cholesky!(A::Hermitian{Float64, Matrix{Float64}}, ::Val{false}; check::Bool)",
" @ LinearAlgebra /Applications/julia-rosetta/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:266",
" [3] cholesky!(A::Matrix{Float64}, ::Val{false}; check::Bool)",
" @ LinearAlgebra /Applications/julia-rosetta/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:298",
" [4] #cholesky#143",
" @ /Applications/julia-rosetta/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:394 [inlined]",
" [5] cholesky (repeats 2 times)",
" @ /Applications/julia-rosetta/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:394 [inlined]",
" [6] top-level scope",
" @ In[14]:1",
" [7] eval",
" @ ./boot.jl:373 [inlined]",
" [8] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
" @ Base ./loading.jl:1196"
]
}
],
"source": [
"cholesky(B) # throws an error"
]
},
{
"cell_type": "markdown",
"id": "4677b31b",
"metadata": {},
"source": [
"It's not hard to manufacture an SPD matrix to try out the Cholesky factorization."
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "7681ce8f",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Cholesky{Float64, Matrix{Float64}}\n",
"U factor:\n",
"4×4 UpperTriangular{Float64, Matrix{Float64}}:\n",
" 11.0905 10.0085 7.57402 8.02486\n",
" ⋅ 6.76973 5.78976 6.30497\n",
" ⋅ ⋅ 6.86388 0.686969\n",
" ⋅ ⋅ ⋅ 2.31886"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"B = A'*A\n",
"cf = cholesky(B)"
]
},
{
"cell_type": "markdown",
"id": "14e39e9d",
"metadata": {},
"source": [
"What's returned is a factorization object. Another step is required to extract the factor as a matrix:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "1323fa0e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 UpperTriangular{Float64, Matrix{Float64}}:\n",
" 11.0905 10.0085 7.57402 8.02486\n",
" ⋅ 6.76973 5.78976 6.30497\n",
" ⋅ ⋅ 6.86388 0.686969\n",
" ⋅ ⋅ ⋅ 2.31886"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"R = cf.U"
]
},
{
"cell_type": "markdown",
"id": "b1f0ca48",
"metadata": {},
"source": [
"Here we validate the factorization:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "57e470d9",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"7.37383175345327e-17"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"opnorm(R'*R - B) / opnorm(B)"
]
},
{
"cell_type": "markdown",
"id": "954ba478",
"metadata": {},
"source": [
"```{raw} html\n",
"

\n",
"```\n",
"\n",
"```{raw} latex\n",
"%%end demo%%\n",
"```\n",
"\n",
"## Exercises\n",
"\n",
"1. ✍ For each matrix, use {eq}`diag-dominant` to determine whether it is diagonally dominant.\n",
"\n",
" ```{math}\n",
" \\mathbf{A} =\n",
" \\begin{bmatrix}\n",
" 3 & 1 & 0 & 1 \\\\\n",
" 0 & -2 & 0 & 1 \\\\\n",
" -1 & 0 & 4 & -1 \\\\\n",
" 0 & 0 & 0 & 6\n",
" \\end{bmatrix},\n",
" \\quad\n",
" \\mathbf{B} =\n",
" \\begin{bmatrix}\n",
" 1 & 0 & 0 & 0 & 0 \\\\\n",
" 0 & 1 & 0 & 0 & 0 \\\\\n",
" 0 & 0 & 1 & 0 & 0 \\\\\n",
" 0 & 0 & 0 & 1 & 0 \\\\\n",
" 0 & 0 & 0 & 0 & 0\n",
" \\end{bmatrix},\n",
" \\quad \\mathbf{C} =\n",
" \\begin{bmatrix}\n",
" 2 & -1 & 0 & 0 \\\\\n",
" -1 & 2 & -1 & 0 \\\\\n",
" 0 & -1 & 2 & -1 \\\\\n",
" 0 & 0 & -1 & 2\n",
" \\end{bmatrix}.\n",
" ```\n",
"\n",
"2. ⌨ For each matrix, use inspection or `cholesky` in Julia to determine whether it is SPD.\n",
"\n",
" ```{math}\n",
" \\mathbf{A} =\n",
" \\begin{bmatrix}\n",
" 1 & 0 & -1 \\\\ 0 & 4 & 5 \\\\ -1 & 5 & 10\n",
" \\end{bmatrix},\n",
" \\qquad\n",
" \\mathbf{B} =\n",
" \\begin{bmatrix}\n",
" 1 & 0 & 1 \\\\ 0 & 4 & 5 \\\\ -1 & 5 & 10\n",
" \\end{bmatrix},\n",
" \\qquad\n",
" \\mathbf{C} =\n",
" \\begin{bmatrix}\n",
" 1 & 0 & 1 \\\\ 0 & 4 & 5 \\\\ 1 & 5 & 1\n",
" \\end{bmatrix}.\n",
" ```\n",
"\n",
"3. ✍ Show that the diagonal entries of a symmetric positive definite matrix are positive numbers. (Hint: Apply certain special cases of {eq}`SPD-def`.)\n",
"\n",
"4. ⌨ Using {numref}`Function {number}