Skip to content

Instantly share code, notes, and snippets.

@tobydriscoll

tobydriscoll/TB-Lecture-23-julia.ipynb

Created October 17, 2016 17:16
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save tobydriscoll/17ce89eb2c827e0b870877219ec64fe8 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save tobydriscoll/17ce89eb2c827e0b870877219ec64fe8 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
TB-Lecture-23-julia.ipynb
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
TB-Lecture-23-matlab.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Lecture 23: Cholesky factorization\n",
"\n",
"## Symmetric elimination\n",
"\n",
"Suppose we have a symmetric (or hermitian) matrix $A$. Let's proceed with Gaussian elimination a little differently."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A =\n",
"\n",
" 33 7 12 17\n",
" 7 23 17 22\n",
" 12 17 13 27\n",
" 17 22 27 3\n"
]
}
],
"source": [
"A = magic(4);\n",
"A = A + A' + eye(4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We start as before, with a unit lower triangular $L_1$ that puts zeros in the first column."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A1 =\n",
"\n",
" 33.0000 7.0000 12.0000 17.0000\n",
" 0 21.5152 14.4545 18.3939\n",
" 0 14.4545 8.6364 20.8182\n",
" 0 18.3939 20.8182 -5.7576\n"
]
}
],
"source": [
"L1 = eye(4); L1(2:4,1) = -A(2:4,1)/A(1,1);\n",
"A1 = L1*A"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"That didn't change the first row of $A$ at all. If we were to transpose $A_1$ and apply $L_1$ on the left again, it would have the same zeroing effect."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n",
"\n",
" 33.0000 0 0 0\n",
" 0 21.5152 14.4545 18.3939\n",
" 0 14.4545 8.6364 20.8182\n",
" 0 18.3939 20.8182 -5.7576\n"
]
}
],
"source": [
"L1*A1'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This result is symmetric, which is no accident. In fact we can get here simply by"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A1 =\n",
"\n",
" 33.0000 0 0 0\n",
" 0 21.5152 14.4545 18.3939\n",
" 0 14.4545 8.6364 20.8182\n",
" 0 18.3939 20.8182 -5.7576\n"
]
}
],
"source": [
"A1 = L1*A*L1'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As long as we apply elimination symmetrically on the left and right, we preserve the symmetry of $A$ as the elimination proceeds. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A2 =\n",
"\n",
" 33.0000 0 0 0\n",
" 0 21.5152 0 0\n",
" 0 0 -1.0746 8.4606\n",
" 0 0 8.4606 -21.4831\n"
]
}
],
"source": [
"L2 = eye(4); L2(3:4,2) = -A1(3:4,2)/A1(2,2);\n",
"A2 = L2*A1*L2'"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A3 =\n",
"\n",
" 33.0000 0 0 0\n",
" 0 21.5152 0 0\n",
" 0 0 -1.0746 -0.0000\n",
" 0 0 0 45.1258\n"
]
}
],
"source": [
"L3 = eye(4); L3(4,3) = -A2(4,3)/A2(3,3);\n",
"A3 = L3*A2*L3'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This matrix is (up to roundoff) symmetric and upper triangular, hence diagonal. In summary, we have a symmetric factorization of $A$:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n",
"\n",
" 2.8422e-14\n"
]
}
],
"source": [
"L = inv(L1)*inv(L2)*inv(L3); D = diag(diag(A3));\n",
"norm( A - L*D*L' )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Like unpivoted LU factorization, the $A=LDL^*$ factorization is unstable. It can be pivoted into stability, but the details are a bit more complex. We're after something more specific, anyhow."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## HPD matrices and the Cholesky factorization\n",
"\n",
"If in addition to symmetry/hermitianness we also know that the quadratic form $x^*Ax$ is positive for all nonzero $x$, then $A$ is said to be **positive definite** (often written **SPD** or **HPD**). By taking $x=e_k$, we can see that $a_{kk}>0$; i.e., all diagonal entries are positive. \n",
"\n",
"Let's rerun the symmetric factorization, with a matrix guaranteed to be HPD."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"B =\n",
"\n",
" 15.7100 9.7000 11.3000 10.9000\n",
" 9.7000 13.5100 12.9000 11.5000\n",
" 11.3000 12.9000 13.3100 10.1000\n",
" 10.9000 11.5000 10.1000 15.1100\n"
]
}
],
"source": [
"B = A'*A/100"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our usual move here is to apply a unit lower triangular matrix. This leaves us with a diagonal term that ended up in the $D$ matrix above. But instead we could distribute that term between the symmetric factors: because $b_{11}>0$, we can define"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"L1 =\n",
"\n",
" 0.2523 0 0 0\n",
" -0.6174 1.0000 0 0\n",
" -0.7193 0 1.0000 0\n",
" -0.6938 0 0 1.0000\n"
]
}
],
"source": [
"c = sqrt(B(1,1));\n",
"L1 = eye(4); \n",
"L1(1,1) = 1/c;\n",
"L1(2:4,1) = -B(2:4,1)/c^2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now the row elimination step still makes the zeros we want..."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n",
"\n",
" 3.9636 2.4473 2.8510 2.7500\n",
" 0 7.5208 5.9229 4.7699\n",
" 0.0000 5.9229 5.1821 2.2598\n",
" 0 4.7699 2.2598 7.5473\n"
]
}
],
"source": [
"L1*B"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"...and the symmetric elimination introduces the zeros _and_ the 1 on the diagonal."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"B1 =\n",
"\n",
" 1.0000 0 0.0000 0\n",
" 0 7.5208 5.9229 4.7699\n",
" 0.0000 5.9229 5.1821 2.2598\n",
" 0 4.7699 2.2598 7.5473\n"
]
}
],
"source": [
"B1 = L1*B*L1'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A few details are between us and an algorithm. First of all, can we continue? The algorithm proceeds recursively (or iteratively, depending on your point of view) if $B_1$ is HPD, so that we have a positive number in the (2,2) pivot spot. This is automatic: if $Z$ is any nonsingular matrix, then\n",
"\n",
"$$ x^*(ZBZ^*)x = (Z^*x)^* B (Z^*x) \\ge 0,$$\n",
"\n",
"and it can be zero only if $Z^*x=0$, which implies $x=0$."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"B2 =\n",
"\n",
" 1.0000 0 0.0000 0\n",
" 0 1.0000 -0.0000 -0.0000\n",
" 0.0000 -0.0000 0.5175 -1.4967\n",
" 0 0 -1.4967 4.5221\n"
]
}
],
"source": [
"c = sqrt(B1(2,2));\n",
"L2 = eye(4); L2(2,2) = 1/c;\n",
"L2(3:4,2) = -B1(3:4,2)/c^2;\n",
"B2 = L2*B1*L2'"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"B3 =\n",
"\n",
" 1.0000 0 0.0000 0.0000\n",
" 0 1.0000 -0.0000 -0.0000\n",
" 0.0000 -0.0000 1.0000 0.0000\n",
" 0.0000 -0.0000 0.0000 0.1939\n",
"\n",
"\n",
"B4 =\n",
"\n",
" 1.0000 0 0.0000 0.0000\n",
" 0 1.0000 -0.0000 -0.0000\n",
" 0.0000 -0.0000 1.0000 0.0000\n",
" 0.0000 -0.0000 0.0000 1.0000\n"
]
}
],
"source": [
"c = sqrt(B2(3,3));\n",
"L3 = eye(4); L3(3,3) = 1/c;\n",
"L3(4:4,3) = -B2(4:4,3)/c^2;\n",
"B3 = L3*B2*L3'\n",
"L4 = eye(4); L4(4,4) = 1/sqrt(B3(4,4));\n",
"B4 = L4*B3*L4'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Another important aspect of the LU factorization is that the terms in $L$ are found easily, because the inverses and products of the row operation matrices are trivial. Now if we write\n",
"\n",
"$$L_1 = \\begin{bmatrix} c^{-1} & 0 \\\\ c^{-2} v & I \\end{bmatrix},$$\n",
"\n",
"then you can quickly check that\n",
"\n",
"$$L_1^{-1} = \\begin{bmatrix} c & 0 \\\\ -c^{-1}v & I \\end{bmatrix}. $$"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n",
"\n",
" 3.9636 0 -0.0000 -0.0000\n",
" 2.4473 1.0000 -0.0000 -0.0000\n",
" 2.8510 0 1.0000 -0.0000\n",
" 2.7500 0 -0.0000 1.0000\n"
]
}
],
"source": [
"inv(L1)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n",
"\n",
" 0.7194\n",
" 13.4833\n",
" 15.7073\n",
" 15.1513\n"
]
}
],
"source": [
"[c; B(2:4,1)/c]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Moreover, thanks to the ordering of the matrices, the nonunit term on the diagonal of one factor only multiplies an identity row to its right. Constructing the complete triangular matrix can therefore be done one column at a time, as before."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n",
"\n",
" 3.9636 0 -0.0000 -0.0000\n",
" 2.4473 1.0000 -0.0000 -0.0000\n",
" 2.8510 0 1.0000 -0.0000\n",
" 2.7500 0 -0.0000 1.0000\n",
"\n",
"\n",
"ans =\n",
"\n",
" 1.0000 0 0 0\n",
" 0 2.7424 0.0000 0.0000\n",
" 0 2.1597 1.0000 0.0000\n",
" 0 1.7393 -0.0000 1.0000\n",
"\n",
"\n",
"ans =\n",
"\n",
" 1.0000 0 0 0\n",
" 0 1.0000 0 0\n",
" 0 0 0.7194 -0.0000\n",
" 0 0 -2.0804 1.0000\n"
]
}
],
"source": [
"inv(L1)\n",
"inv(L2)\n",
"inv(L3)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n",
"\n",
" 3.9636 -0.0000 0.0000 -0.0000\n",
" 2.4473 2.7424 0.0000 -0.0000\n",
" 2.8510 2.1597 0.7194 0.0000\n",
" 2.7500 1.7393 -2.0804 0.4403\n"
]
}
],
"source": [
"inv(L1)*inv(L2)*inv(L3)*inv(L4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The resulting **Cholesky factorization** takes the form $A=R^*R$ for an upper triangular $R$ with positive diagonal entries. Of course we don't need to compute both an upper and a lower triangle, and as a result it takes asympototically half as many flops as standard $LU$ factorization. Furthermore, the factorization is backward stable _without_ pivoting."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Matlab",
"language": "matlab",
"name": "matlab"
},
"language_info": {
"codemirror_mode": "octave",
"file_extension": ".m",
"help_links": [
{
"text": "MetaKernel Magics",
"url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
}
],
"mimetype": "text/x-octave",
"name": "matlab",
"version": "0.11.0"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment