Skip to content

Instantly share code, notes, and snippets.

@csullivan

csullivan/qr_decomp.ipynb

Created October 1, 2017 22:06
Show Gist options
  • Save csullivan/dccb66b2bb54796c14142f619d2b1761 to your computer and use it in GitHub Desktop.
Save csullivan/dccb66b2bb54796c14142f619d2b1761 to your computer and use it in GitHub Desktop.
Download ZIP
Matrix QR Decompositions
Display the source blob
Display the rendered blob
Raw
qr_decomp.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# QR decompositions - Chris Sullivan - 10/02/17"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This iPython notebook compares the QR decomposition of two matrices using 4 different methods: \n",
" \n",
"* Classical Gram-Schmidt Orthogonalization \n",
"* Modified Gram-Schmidt Orthogonalization \n",
"* Householder Orthogonalization \n",
"* numpy.linalg.qr - NumPy's built in QR decomposition (householder) \n",
" \n",
"In both case 1 and case 2, we see that the Gram-Schmidt and Householder methods differ by a minus sign, which suggest the solution is unique up to the sign. In comparing the orthogonality of the resulting unitary matrix Q for case 2, we find that both of the Gram-Schmidt techniques have more error and are thus less orthogonal than in the case of the Householder orthogonalization. In comparing the performance, the classical Gram-Schmidt executes the fastest, of the python-implemented routines. The numpy householder implementation utilizes compiled LAPACK libraries and is thus faster than each of the implementations in this work. However, I will note that the Householder implementation for this assignment\n",
"produced the best orthogonality out of all the methods, including the numpy implementation, this could imply that the numpy implementation is sacrificing a small amount of precision for performance. "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import numpy as np\n",
"import numpy.linalg"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"A1 = np.array([ # case 1\n",
" [1,2,3],\n",
" [4,5,6],\n",
" [7,8,7],\n",
" [4,2,3],\n",
" [4,2,2]\n",
"], dtype=np.float64)\n",
"\n",
"A2 = np.array([ # case 2\n",
" [0.70000,0.70711],\n",
" [0.70001,0.70711]\n",
"], dtype=np.float64)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Classic Gram-Schmidt QR decomposition"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def classic_gram_schmidt(A):\n",
" \"\"\" The jth column of Q is the jth column \n",
" of A with the projected subtractions of \n",
" every preceeding column in A, then normalized \"\"\"\n",
" \n",
" nrows,ncols = A.shape\n",
" \n",
" Q = np.copy(A,order='F') # columns in contiguous memory\n",
" R = np.zeros(shape=(ncols,ncols))\n",
" for colj in range(ncols):\n",
" for coli in range(colj):\n",
" R[coli,colj] = Q[:,coli].dot(Q[:,colj])\n",
" Q[:,colj] -= R[coli,colj]*Q[:,coli]\n",
" norm2 = np.linalg.norm(Q[:,colj],ord=2)\n",
" if norm2:\n",
" Q[:,colj] /= norm2\n",
" else:\n",
" raise RuntimeError(\"One of the vectors of A is the 0 vector, or A is not of full rank\")\n",
" \n",
" R[colj,colj] = norm2\n",
" \n",
" return Q,R"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Modified Gram-Schmidt QR decomposition"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def modified_gram_schmidt(A):\n",
" \"\"\" Subtract the ith normalized column of A from every jth column of A such that j>i \"\"\"\n",
" nrows,ncols = A.shape\n",
" \n",
" Q = np.copy(A,order='F') # columns in contiguous memory\n",
" R = np.zeros(shape=(ncols,ncols))\n",
" \n",
" for coli in range(ncols):\n",
" R[coli,coli] = np.linalg.norm(Q[:,coli],ord=2)\n",
" if R[coli,coli]:\n",
" Q[:,coli] /= R[coli,coli]\n",
" else: \n",
" raise RuntimeError(\"One of the columns of A is the 0 vector\")\n",
" for colj in range(coli+1,ncols):\n",
" R[coli,colj] = Q[:,coli].dot(Q[:,colj])\n",
" Q[:,colj] -= R[coli,colj]*Q[:,coli] \n",
" return Q,R"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Householder algorithm for QR decomposition"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def householder_reflector(x):\n",
" v = np.copy(x)\n",
" v[0] += np.sign(x[0])*np.linalg.norm(x,ord=2)\n",
" return v\n",
"\n",
"def householder_qr(A):\n",
" R = np.copy(A,order='F')\n",
" nrows,ncols = R.shape\n",
" \n",
" hvectors = []\n",
" for icol in range(ncols):\n",
" v = householder_reflector(R[icol:,icol])\n",
" norm2 = np.linalg.norm(v,ord=2)\n",
" if norm2:\n",
" v /= norm2\n",
" else:\n",
" raise RuntimeError(\"Norm of zero vector\")\n",
" \n",
" R[icol:,icol:] -= 2*np.outer(v,v.dot(R[icol:,icol:]))\n",
" hvectors.append(v)\n",
" \n",
" # calculate Q according to 10.3\n",
" Q = np.eye(nrows)\n",
" for icol in range(ncols):\n",
" for k in reversed(range(ncols)): \n",
" Q[k:,icol] -= 2*hvectors[k]*(hvectors[k].dot(Q[k:,icol]))\n",
"\n",
" return Q[:nrows,:ncols], R[:ncols,:ncols]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Problem 2: Case I"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 0.10101525 0.31617307 0.5419969 ]\n",
" [ 0.40406102 0.3533699 0.51618752]\n",
" [ 0.70710678 0.39056673 -0.52479065]\n",
" [ 0.40406102 -0.55795248 0.38714064]\n",
" [ 0.40406102 -0.55795248 -0.12044376]]\n",
"\n",
"[[ 9.89949494 9.49543392 9.69746443]\n",
" [ 0. 3.29191961 3.01294337]\n",
" [ 0. 0. 1.97011572]]\n",
"\n",
"CPU times: user 0 ns, sys: 0 ns, total: 0 ns\n",
"Wall time: 250 µs\n",
"\n",
"\n",
"[[ 0.10101525 0.31617307 0.5419969 ]\n",
" [ 0.40406102 0.3533699 0.51618752]\n",
" [ 0.70710678 0.39056673 -0.52479065]\n",
" [ 0.40406102 -0.55795248 0.38714064]\n",
" [ 0.40406102 -0.55795248 -0.12044376]]\n",
"\n",
"[[ 9.89949494 9.49543392 9.69746443]\n",
" [ 0. 3.29191961 3.01294337]\n",
" [ 0. 0. 1.97011572]]\n",
"\n",
"CPU times: user 0 ns, sys: 0 ns, total: 0 ns\n",
"Wall time: 248 µs\n",
"\n",
"\n",
"[[-0.10101525 -0.31617307 0.5419969 ]\n",
" [-0.40406102 -0.3533699 0.51618752]\n",
" [-0.70710678 -0.39056673 -0.52479065]\n",
" [-0.40406102 0.55795248 0.38714064]\n",
" [-0.40406102 0.55795248 -0.12044376]]\n",
"\n",
"[[ -9.89949494e+00 -9.49543392e+00 -9.69746443e+00]\n",
" [ -8.88178420e-16 -3.29191961e+00 -3.01294337e+00]\n",
" [ -8.88178420e-16 0.00000000e+00 1.97011572e+00]]\n",
"\n",
"CPU times: user 0 ns, sys: 0 ns, total: 0 ns\n",
"Wall time: 744 µs\n",
"\n",
"\n",
"[[-0.10101525 -0.31617307 0.5419969 ]\n",
" [-0.40406102 -0.3533699 0.51618752]\n",
" [-0.70710678 -0.39056673 -0.52479065]\n",
" [-0.40406102 0.55795248 0.38714064]\n",
" [-0.40406102 0.55795248 -0.12044376]]\n",
"\n",
"[[-9.89949494 -9.49543392 -9.69746443]\n",
" [ 0. -3.29191961 -3.01294337]\n",
" [ 0. 0. 1.97011572]]\n",
"\n",
"CPU times: user 0 ns, sys: 0 ns, total: 0 ns\n",
"Wall time: 191 µs\n",
"\n",
"\n"
]
}
],
"source": [
"for qr_decomp in [classic_gram_schmidt, modified_gram_schmidt, householder_qr, np.linalg.qr]:\n",
" Q,R = qr_decomp(A1)\n",
" print(Q)\n",
" print()\n",
" print(R)\n",
" print()\n",
" %time qr_decomp(A1)\n",
" print('\\n')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Problem 2: Case II"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 0.70710173 0.70711183]\n",
" [ 0.70711183 -0.70710173]]\n",
"\n",
"[[ 9.89956565e-01 1.00000455e+00]\n",
" [ 0.00000000e+00 7.14283864e-06]]\n",
"\n",
"Orthogonality check:\n",
"[[ 0.00000000e+00 2.30143682e-11]\n",
" [ 2.30143682e-11 0.00000000e+00]]\n",
"\n",
"\n",
"[[ 0.70710173 0.70711183]\n",
" [ 0.70711183 -0.70710173]]\n",
"\n",
"[[ 9.89956565e-01 1.00000455e+00]\n",
" [ 0.00000000e+00 7.14283864e-06]]\n",
"\n",
"Orthogonality check:\n",
"[[ 0.00000000e+00 2.30143682e-11]\n",
" [ 2.30143682e-11 0.00000000e+00]]\n",
"\n",
"\n",
"[[-0.70710173 0.70711183]\n",
" [-0.70711183 -0.70710173]]\n",
"\n",
"[[ -9.89956565e-01 -1.00000455e+00]\n",
" [ 0.00000000e+00 7.14283864e-06]]\n",
"\n",
"Orthogonality check:\n",
"[[ 0. 0.]\n",
" [ 0. 0.]]\n",
"\n",
"\n",
"[[-0.70710173 -0.70711183]\n",
" [-0.70711183 0.70710173]]\n",
"\n",
"[[ -9.89956565e-01 -1.00000455e+00]\n",
" [ 0.00000000e+00 -7.14283864e-06]]\n",
"\n",
"Orthogonality check:\n",
"[[ -2.22044605e-16 -5.55111512e-17]\n",
" [ -5.55111512e-17 0.00000000e+00]]\n",
"\n",
"\n"
]
}
],
"source": [
"for qr_decomp in [classic_gram_schmidt, modified_gram_schmidt, householder_qr, np.linalg.qr]:\n",
" Q,R = qr_decomp(A2)\n",
" print(Q)\n",
" print()\n",
" print(R)\n",
" print()\n",
" print(\"Orthogonality check:\")\n",
" print(Q.T.dot(Q) - np.eye(Q.shape[0]))\n",
" print('\\n')\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The householder algorithm has the best orthogonality, even better than the built in numpy QR decomposition."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.4.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment