\n",
"```\n",
"\n",
"```{raw} latex\n",
"%%start demo%%\n",
"```\n",
"\n",
"```{index} ! Julia; cond\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",
"Julia has a function `cond` to compute matrix condition numbers. By default, the 2-norm is used. As an example, the family of *Hilbert matrices* is famously badly conditioned. Here is the $6\\times 6$ case.\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",
"Type `\\kappa` followed by `Tab` to get the Greek letter $\\kappa$.\n",
"```{raw} latex\n",
"\\end{mdframed}\\end{minipage}\n",
"```\n",
"::::"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "1b2c28a2",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5.109816296610315e7"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = [ 1/(i+j) for i in 1:6, j in 1:6 ]\n",
"κ = cond(A)"
]
},
{
"cell_type": "markdown",
"id": "cab75e94",
"metadata": {},
"source": [
"Because $\\kappa\\approx 10^8$, it's possible to lose nearly 8 digits of accuracy in the process of passing from $\\mathbf{A}$ and $\\mathbf{b}$ to $\\mathbf{x}$. That fact is independent of the algorithm; it's inevitable once the data are expressed in finite precision. \n",
"\n",
"Let's engineer a linear system problem to observe the effect of a perturbation. We will make sure we know the exact answer."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "6f1548c3",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6-element Vector{Float64}:\n",
" 4.407142857142857\n",
" 3.564285714285714\n",
" 3.013095238095238\n",
" 2.6174603174603175\n",
" 2.317279942279942\n",
" 2.0807359307359308"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = 1:6\n",
"b = A*x"
]
},
{
"cell_type": "markdown",
"id": "a6092eb6",
"metadata": {},
"source": [
"Now we perturb the system matrix and vector randomly by $10^{-10}$ in norm."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "c039295a",
"metadata": {},
"outputs": [],
"source": [
"ΔA = randn(size(A)); ΔA = 1e-10*(ΔA/opnorm(ΔA));\n",
"Δb = randn(size(b)); Δb = 1e-10*normalize(Δb);"
]
},
{
"cell_type": "markdown",
"id": "668775d4",
"metadata": {},
"source": [
"We solve the perturbed problem using pivoted LU and see how the solution was changed."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "786772d3",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6-element Vector{Float64}:\n",
" -5.471659263400763e-5\n",
" 0.0009283239117472419\n",
" -0.004843432603480302\n",
" 0.010718598499813048\n",
" -0.010580070933729147\n",
" 0.0038388020192696715"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"new_x = ((A+ΔA) \\ (b+Δb))\n",
"Δx = new_x - x"
]
},
{
"cell_type": "markdown",
"id": "af412cd8",
"metadata": {},
"source": [
"Here is the relative error in the solution."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "7549d109",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"relative_error = norm(Δx) / norm(x) = 0.0017093353064577405\n"
]
}
],
"source": [
"@show relative_error = norm(Δx) / norm(x);"
]
},
{
"cell_type": "markdown",
"id": "fec1226d",
"metadata": {},
"source": [
"And here are upper bounds predicted using the condition number of the original matrix."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "4cd273ec",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Upper bound due to b: 0.0006723667712613617\n",
"Upper bound due to A: 0.004566989872745153\n"
]
}
],
"source": [
"println(\"Upper bound due to b: $(κ*norm(Δb)/norm(b))\")\n",
"println(\"Upper bound due to A: $(κ*opnorm(ΔA)/opnorm(A))\")"
]
},
{
"cell_type": "markdown",
"id": "e4030fa2",
"metadata": {},
"source": [
"Even if we didn't make any manual perturbations to the data, machine roundoff does so at the relative level of $\\macheps$."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "94073540",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"relative_error = norm(Δx) / norm(x) = 1.2278255878021855e-9\n",
"rounding_bound = κ * eps() = 1.134607140820324e-8\n"
]
}
],
"source": [
"Δx = A\\b - x\n",
"@show relative_error = norm(Δx) / norm(x);\n",
"@show rounding_bound = κ*eps();"
]
},
{
"cell_type": "markdown",
"id": "17ac847c",
"metadata": {},
"source": [
"Larger Hilbert matrices are even more poorly conditioned:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "ece92a79",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3.870898056087325e17"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = [ 1/(i+j) for i=1:14, j=1:14 ];\n",
"κ = cond(A)"
]
},
{
"cell_type": "markdown",
"id": "fe586c4f",
"metadata": {},
"source": [
"Note that $\\kappa$ exceeds $1/\\macheps$. In principle we therefore may end up with an answer that has relative error greater than 100%."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "5b9aa3a8",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"85.95120295689817"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rounding_bound = κ*eps()"
]
},
{
"cell_type": "markdown",
"id": "d797b9d3",
"metadata": {},
"source": [
"Let's put that prediction to the test."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "66521b4a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"relative_error = norm(Δx) / norm(x) = 2.3271868502658917\n"
]
}
],
"source": [
"x = 1:14\n",
"b = A*x \n",
"Δx = A\\b - x\n",
"@show relative_error = norm(Δx) / norm(x);"
]
},
{
"cell_type": "markdown",
"id": "47232167",
"metadata": {},
"source": [
"As anticipated, the solution has zero accurate digits in the 2-norm.\n",
"```{raw} html\n",
"

\n",
"```\n",
"\n",
"```{raw} latex\n",
"%%end demo%%\n",
"```\n",
"\n",
"## Residual and backward error\n",
"\n",
"Suppose that $\\mathbf{A}\\mathbf{x}=\\mathbf{b}$ and $\\tilde{\\mathbf{x}}$ is a computed estimate of the solution $\\mathbf{x}$. The most natural quantity to study is the error, $\\mathbf{x}-\\tilde{\\mathbf{x}}$. Normally we can't compute it because we don't know the exact solution. However, we can compute something related.\n",
"\n",
"```{index} ! residual; of a linear system\n",
"```\n",
"\n",
"::::{proof:definition} Residual of a linear system\n",
"For the problem $\\mathbf{A}\\mathbf{x}=\\mathbf{b}$, the **residual** at a solution estimate $\\tilde{\\mathbf{x}}$ is \n",
"\n",
"```{math}\n",
" :label: residual\n",
" \\mathbf{r} = \\mathbf{b} - \\mathbf{A}\\tilde{\\mathbf{x}}.\n",
"```\n",
"::::\n",
"\n",
"```{index} backward error; in a linear system\n",
"```\n",
"\n",
"Obviously, a zero residual means that $\\tilde{\\mathbf{x}}=\\mathbf{x}$, and we have the exact solution. What happens more generally? Note that $\\mathbf{A}\\tilde{\\mathbf{x}}=\\mathbf{b}-\\mathbf{r}$. That is, $\\tilde{\\mathbf{x}}$ solves the linear system problem for a right-hand side that is changed by $-\\mathbf{r}$. This is precisely what is meant by backward error.\n",
"\n",
"Hence residual and backward error are the same thing for a linear system. What is the connection to the (forward) error? We can reconnect with {eq}`linsyscondb` by the definition $\\mathbf{h} = \\tilde{\\mathbf{x}}-\\mathbf{x}$, in which case \n",
"\n",
"$$\\mathbf{d} = \\mathbf{A}(\\mathbf{x}+\\mathbf{h})-\\mathbf{b}=\\mathbf{A}\\mathbf{h} = -\\mathbf{r}.$$ \n",
"\n",
"Thus {eq}`linsyscondb` is equivalent to\n",
"\n",
"```{math}\n",
" :label: residualcond\n",
" \\frac{\\| \\mathbf{x}-\\tilde{\\mathbf{x}} \\|}{\\| \\mathbf{x} \\|} \\le\n",
" \\kappa(\\mathbf{A}) \\frac{\\| \\mathbf{r} \\|}{\\| \\mathbf{b} \\|}.\n",
"```\n",
"\n",
"Equation {eq}`residualcond` says that the gap between relative error and the relative residual is a multiplication by the matrix condition number. \n",
"\n",
"```{proof:observation}\n",
"When solving a linear system, all that can be expected is that the backward error, not the error, is small.\n",
"```\n",
"\n",
"## Exercises\n",
"\n",
"1. ⌨ Refer to {numref}`Demo {number}