{
"metadata": {
"language": "Julia",
"name": ""
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Julia Issue Note #5705\n",
"## Linear algebra test failures for linear solve with triangular matrices\n",
"\n",
"Recent Julia issues like [#5705](https://github.com/JuliaLang/julia/pull/5705) have adjusted the bounds on linear algebra tests. The purpose of this notebook is to provide more extensive numerical demonstrations of why the current linear algebra tests are fundamentally broken and need to be systematically replaced by more theoretically justified error bounds, as explained in [#5605](https://github.com/JuliaLang/julia/issues/5605).\n",
"\n",
"## The current testing strategy\n",
"\n",
"Here is a snippet of code representing a typical [linear algebra test](https://github.com/JuliaLang/julia/blob/master/test/linalg.jl) in Julia:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#Test linear solve on upper triangular matrices\n",
"n=10\n",
"eltya=Float64 #eltype of matrix\n",
"eltyb=Float64 #eltype of vector\n",
"\n",
"A=convert(Matrix{eltya}, randn(n,n))\n",
"A=triu(A)\n",
"b=convert(Vector{eltyb}, randn(n))\n",
"\n",
"\u03b5 = max(eps(abs(float(one(eltya)))),eps(abs(float(one(eltyb)))))\n",
"x = A \\ b\n",
"#@test_approx_eq_eps triu(a)*x b 20000\u03b5\n",
"#Instead of the actual assertion, show the values being compared\n",
"@show norm(A*x), norm(b), 20000\u03b5\n",
"abs(norm(A*x) - norm(b))/\u03b5*1"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"(norm(*(A,x)),norm(b),*(20000,\u03b5)) => (4.026090098980171,4.026090098980171,4.440892098500626e-12)"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n"
]
},
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 1,
"text": [
"0.0"
]
}
],
"prompt_number": 1
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Mathematically, the test being performed corresponds to asserting that the inequality\n",
"\n",
"$$\n",
"\\left\\vert \\left\\Vert A x\\right\\Vert - \\left\\Vert b \\right\\Vert \\right\\vert \\le C \\epsilon\n",
"$$\n",
"\n",
"is valid, where $\\epsilon$ is machine epsilon and $C = 20000$ is\n",
"an arbitrary magic constant.\n",
"\n",
"Most of the time, this test will appear to be well-behaved. The deviation between $Ax$ and $b$ appears to be small in norm, and it looks like the magic number 20000 would _obviously_ cover any contingency. But how robust is this observation?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Distribution of magic numbers"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"t=10^6\n",
"n=10\n",
"T=Float64\n",
"\n",
"b=randn(n)\n",
"c=zeros(t) #log10 of magic numbers\n",
"\u03b5=eps(real(one(T)))\n",
"for i=1:t\n",
" A = convert(Matrix{T}, triu(randn(n,n)))\n",
" x = A \\ b\n",
" c[i] = max(log10(abs(norm(A*x)-norm(b))/\u03b5), -2) #arbitrary floor for plotting 0 on log scale\n",
"end"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"@show minimum(c), maximum(c)\n",
"x, y = hist(c, 50)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"(minimum(c),maximum(c)) => "
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"(-2.0,12.659850021280592)\n"
]
},
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 3,
"text": [
"(-2.5:0.5:13.0,[67288,0,0,0,0,99814,180104,156724,160350,120193 \u2026 65,24,14,0,5,0,0,0,0,1])"
]
}
],
"prompt_number": 3
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The statistics show that every now and then, a nearly singular matrix is sampled, causing the magic number to blow up.\n",
"\n",
"Where do the old and new bounds $15000\\epsilon$ and $20000\\epsilon$ lie on the distribution?"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"using Gadfly\n",
"plot(x=x, y=y/(t*(x[2]-x[1])), Geom.bar, Guide.XLabel(\"log10(c)\"), Guide.YLabel(\"Density\"), \n",
" xintercept=[log10(15000), log10(20000)], Geom.vline(color=\"orange\"))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
""
],
"metadata": {},
"output_type": "display_data"
},
{
"html": [
""
],
"metadata": {},
"output_type": "display_data"
},
{
"html": [
"\n",
"\n"
],
"metadata": {},
"output_type": "display_data",
"text": [
"D3(120.0,80.0,IOBuffer([0x66,0x75,0x6e,0x63,0x74,0x69,0x6f,0x6e,0x20,0x64 \u2026 0x20,0x3d,0x20,0x64,0x72,0x61,0x77,0x0a,0x7d,0x0a],true,true,true,false,26051,9223372036854775807,26052),0,[],[],0,3,Dict{Uint64,(Any,Int64)}(),true,false,nothing,true)"
]
},
{
"html": [],
"metadata": {},
"output_type": "pyout",
"prompt_number": 4,
"text": [
"Plot(...)"
]
}
],
"prompt_number": 4
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Despite the large change in the absolute magnitude of the magic constant, we see that this only shifts the threshold slightly on this histogram of $log(c)$, thanks to the extremely long tail in the magic number $c$. What is the probability that a given random matrix will fail the old and new cutoffs?"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"for cutoff in (15000, 20000)\n",
" p = sum([c .> log10(cutoff)])/t\n",
" @show cutoff, p\n",
"end"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"(cutoff,p) => "
]
},
{
"output_type": "stream",
"stream": "stderr",
"text": [
"\n",
"WARNING: deprecated syntax \"x[i:]\" at /Users/jiahao/.julia/v0.3/DataFrames/src/formula.jl:58.\n",
"Use \"x[i:end]\" instead.\n",
"\n",
"WARNING: deprecated syntax \"x[i:]\" at /Users/jiahao/.julia/v0.3/DataFrames/src/formula.jl:71.\n",
"Use \"x[i:end]\" instead.\n",
"\n",
"WARNING: deprecated syntax \"x[i:]\" at /Users/jiahao/.julia/v0.3/DataFrames/src/formula.jl:77.\n",
"Use \"x[i:end]\" instead.\n",
"\n",
"WARNING: deprecated syntax \"x[i:]\" at /Users/jiahao/.julia/v0.3/DataFrames/src/formula.jl:92.\n",
"Use \"x[i:end]\" instead.\n",
"\n",
"WARNING: deprecated syntax \"x[i:]\" at /Users/jiahao/.julia/v0.3/DataFrames/src/formula.jl:106.\n",
"Use \"x[i:end]\" instead.\n",
"\n",
"WARNING: deprecated syntax \"x[i:]\" at /Users/jiahao/.julia/v0.3/DataFrames/src/formula.jl:111.\n",
"Use \"x[i:end]\" instead.\n",
"\n",
"WARNING: deprecated syntax \"x[i:]\" at /Users/jiahao/.julia/v0.3/DataFrames/src/formula.jl:122.\n",
"Use \"x[i:end]\" instead.\n"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"(15000,0.032434)\n",
"(cutoff,p) => (20000,0.027863)\n"
]
}
],
"prompt_number": 5
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The change of the magic number in PR #5707 decreases the probability of failure by an essentially negligible amount.\n",
"\n",
"Plotting the failure rate vs the magic constant $log(c)$ shows that in the grand scheme of things, the change in the magic constant does very little, and it would take an enormous shift in $c$ to make the failure rate negligible, so much so that the test would be practically useless."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"plot(x=x, y=1-cumsum(y)/t, Geom.line, Guide.XLabel(\"log10(c)\"), Guide.YLabel(\"Probability of failure\"),\n",
"xintercept=[log10(15000), log10(20000)], Geom.vline(color=\"Orange\"))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"\n",
"\n"
],
"metadata": {},
"output_type": "display_data",
"text": [
"D3(120.0,80.0,IOBuffer([0x66,0x75,0x6e,0x63,0x74,0x69,0x6f,0x6e,0x20,0x64 \u2026 0x20,0x3d,0x20,0x64,0x72,0x61,0x77,0x0a,0x7d,0x0a],true,true,true,false,16365,9223372036854775807,16366),0,[],[],0,3,Dict{Uint64,(Any,Int64)}(),true,false,nothing,true)"
]
},
{
"html": [],
"metadata": {},
"output_type": "pyout",
"prompt_number": 6,
"text": [
"Plot(...)"
]
}
],
"prompt_number": 6
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Theoretical error bounds on linear solvers\n",
"\n",
"Let's have a look at which method is being dispatched by the linear solver."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"@which A\\b\n",
"@which \\(Triangular(A), b)\n",
"@which Base.LinAlg.naivesub!(Triangular(A), b)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"naivesub!(A::Triangular{T<:Number},b::AbstractArray{T,1}) at linalg/triangular.jl:127"
],
"metadata": {},
"output_type": "pyout",
"prompt_number": 7,
"text": [
"naivesub!(A::Triangular{T<:Number},b::AbstractArray{T,1}) at linalg/triangular.jl:127"
]
}
],
"prompt_number": 7
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For a triangular matrix $A$, this either dispatches to the LAPACK [xTRTRS](http://www.netlib.org/lapack/explore-html/d6/d6f/dtrtrs_8f.html) subroutine, which in turn calls the substitution solver [xTRSM](http://www.netlib.org/lapack/explore-html/de/da7/dtrsm_8f.html), or the na\u00efve substitution solver `Base.LinAlg.naivesub!`, depending on the element type.\n",
"\n",
"What does theory tell us about the error bounds on substitution? You can derive an error bound using first-order matrix perturbation theory, or if you're lazy, consult a reference like Higham's [_Accuracy and Stability of Numerical Algorithms_](http://epubs.siam.org/doi/book/10.1137/1.9780898718027), Ch. 8 [(pdf)](http://epubs.siam.org/doi/pdf/10.1137/1.9780898718027.ch8)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Error bound results usually come in two flavors. First, we have forward error bounds, which bound the norm of the difference between the computed solution $\\hat{x}$ and the exact solution $x$.\n",
"\n",
"For the triangular system, we have the forward error bound (Higham, Theorem 8.5):\n",
"\n",
"$$\n",
"\\frac{\\left\\Vert x - \\hat{x} \\right\\Vert_\\infty}{\\left\\Vert x \\right\\Vert_\\infty} \\le \\frac{\\textrm{cond}(A,x) \\gamma_n}{1 - \\textrm{cond}(T) \\gamma_n}\n",
"$$\n",
"\n",
"with constant\n",
"\n",
"$$\n",
"\\gamma_n = \\frac{n\\epsilon}{1-n\\epsilon}\n",
"$$\n",
"\n",
"and the Skeel condition numbers\n",
"\n",
"$$\n",
"\\textrm{cond}(A,x) = \\frac{\\left\\Vert \\left\\vert A \\right\\vert \\left\\vert A^{-1} \\right\\vert \\left\\vert x \\right\\vert \\right\\Vert_\\infty}{\\left\\Vert x \\right\\Vert_\\infty}\n",
"$$\n",
"\n",
"and\n",
"\n",
"$$\n",
"\\textrm{cond}(A) = \\left\\Vert \\left\\vert A \\right\\vert \\left\\vert A^{-1} \\right\\vert \\right\\Vert_\\infty\n",
"$$\n",
"\n",
"where $|A|$ means that the modulus is taken componentwise for every element of $A$.\n",
"\n",
"Second, we also have backward error bounds, which assert that the computed solution $\\hat{x} = A\\backslash b$ is an exact solution to a slightly perturbed problem $(A + \\delta A) \\hat{x} = b + \\delta b$. Sometimes results are also known if only $A$ or only $b$ is assumed to be perturbed.\n",
"\n",
"For the triangular system, we have the backward error bound (Higham, Theorem 8.3)\n",
"\n",
"$$\n",
"\\left\\vert \\delta A \\right\\vert \\le \\gamma_n \\left\\vert A \\right\\vert\n",
"$$\n",
"\n",
"where only the matrix $A$ is perturbed."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Forward error bound\n",
"\n",
"Ordinarily, this error bound is useless in practice since the exact solution $x$ is not known. However, thanks to Julia's `BigFloat`s, we can _approximate_ the exact solution with that returned by `big(A) \\ big(b)`."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#Define the Skeel condition numbers\n",
"cond_skeel(A::AbstractMatrix) = norm(abs(inv(A))*abs(A), Inf)\n",
"cond_skeel(A::AbstractMatrix, x::AbstractVector) = norm(abs(inv(A))*abs(A)*abs(x), Inf)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 8,
"text": [
"cond_skeel (generic function with 2 methods)"
]
}
],
"prompt_number": 8
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"t=10^4\n",
"n=10\n",
"T=Float64\n",
"\n",
"b=randn(n)\n",
"c=zeros(t) #log10 of magic numbers\n",
"r=zeros(t) #ratio of lhs/rhs in inequality\n",
"\u03b5=eps(real(one(T)))\n",
"\u03b3 = n*\u03b5/(1-n*\u03b5) #\u03b3_n\n",
"for i=1:t\n",
" A = convert(Matrix{T}, triu(randn(n,n)))\n",
" \u0302x = A \\ b #we called this x earlier\n",
" c[i] = max(log10(abs(norm(A* \u0302x)-norm(b))/\u03b5), -2) #arbitrary floor for plotting 0 on log scale\n",
" bigA = big(A)\n",
" x = bigA \\ b #x is now the exact solution\n",
" lhs = norm(\u0302x-x,Inf)/norm(x,Inf) \n",
" rhs = cond_skeel(bigA, x)*\u03b3/(1-cond_skeel(bigA)*\u03b3)\n",
" r[i] = lhs/rhs\n",
"end"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 9
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"plot(x=c, y=r, Guide.XLabel(\"log10(c)\"), Guide.YLabel(\"lhs/rhs\"))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"\n",
"\n"
],
"metadata": {},
"output_type": "display_data",
"text": [
"D3(120.0,80.0,IOBuffer([0x66,0x75,0x6e,0x63,0x74,0x69,0x6f,0x6e,0x20,0x64 \u2026 0x20,0x3d,0x20,0x64,0x72,0x61,0x77,0x0a,0x7d,0x0a],true,true,true,false,376559,9223372036854775807,376560),0,[],[],1,3,[0x168a35466291c707=>([62.30299170957017,70.56540663720602,60.69215937679767,62.30299170957017,66.9855050688087,65.79837706158487,71.07066107813712,48.37918867924529,83.72324855565678,69.58969156768646 \u2026 62.30299170957017,83.51690222754041,67.87555601388665,76.83651350096237,70.39318416987234,48.37918867924529,66.46693288473628,64.43239622641511,48.37918867924529,68.81783252350775],0),0x941bf19eeee15b5a=>([56.39069659997892,57.89769100842021,57.636525181762146,57.72706088127445,57.830449363170914,57.89735391834255,57.87304390573219,57.86976154390013,57.89924060452837,57.843023136107 \u2026 57.487544433598586,57.899131627069515,57.86191083934261,57.897928771745576,57.881742868818094,51.768343231987515,57.884055517209674,57.891578498215125,55.370008205094926,57.885605998218935],1)],true,false,nothing,true)"
]
},
{
"html": [],
"metadata": {},
"output_type": "pyout",
"prompt_number": 10,
"text": [
"Plot(...)"
]
}
],
"prompt_number": 10
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There doesn't seem to be much correlation between the originally coded test and the test using the forward error bound. But at least the forward error bound test seems to be robust; every point has `lhs <= rhs`."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"x, y=hist(r, 100)\n",
"plot(x=x, y=y/(t*(x[2]-x[1])), Geom.bar, Guide.XLabel(\"lhs/rhs\"), Guide.YLabel(\"Density\"),\n",
"xintercept=[1.0], Geom.vline(color=\"orange\"))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"\n",
"\n"
],
"metadata": {},
"output_type": "display_data",
"text": [
"D3(120.0,80.0,IOBuffer([0x66,0x75,0x6e,0x63,0x74,0x69,0x6f,0x6e,0x20,0x64 \u2026 0x20,0x3d,0x20,0x64,0x72,0x61,0x77,0x0a,0x7d,0x0a],true,true,true,false,44478,9223372036854775807,44479),0,[],[],0,3,Dict{Uint64,(Any,Int64)}(),true,false,nothing,true)"
]
},
{
"html": [],
"metadata": {},
"output_type": "pyout",
"prompt_number": 11,
"text": [
"Plot(...)"
]
}
],
"prompt_number": 11
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"mean(r), std(r)/sqrt(t-1) #mean with its standard error"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 12,
"text": [
"(8.240937806544264e-5,3.067119378454114e-6)"
]
}
],
"prompt_number": 12
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The inequality seems to be fairly loose, but at least nothing violates it."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Backward error bound\n",
"\n",
"From the implicit definition of the perturbed matrix $\\delta A$ in\n",
"\n",
"$$\n",
"(A + \\delta A) \\hat{x} = b\n",
"$$\n",
"\n",
"we can rearrange this to get\n",
"\n",
"$$\n",
"\\delta A \\hat{x} = b - A \\hat{x}\n",
"$$\n",
"\n",
"And taking the norm gives an inequality on the residual\n",
"$$\n",
"\\left\\Vert b - A \\hat{x} \\right\\Vert\n",
"= \\left\\Vert \\delta A \\hat{x} \\right\\Vert\n",
"= \\left\\Vert \\left\\vert \\delta A \\hat{x} \\right\\vert \\right\\Vert\n",
"\\le \\gamma_n \\left\\Vert \\left\\vert A \\right\\vert \\left\\vert \\hat{x} \\right\\vert \\right\\Vert\n",
"\\le \\gamma_n \\left\\Vert A \\right\\Vert \\left\\Vert \\hat{x} \\right\\Vert.\n",
"$$"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"t=10^4\n",
"n=10\n",
"T=Float64\n",
"\n",
"b=randn(n)\n",
"c=zeros(t) #log10 of magic numbers\n",
"r=zeros(t) #ratio of lhs/rhs in inequality\n",
"\u03b5=eps(real(one(T)))\n",
"\u03b3 = n*\u03b5/(1-n*\u03b5) #\u03b3_n\n",
"for i=1:t\n",
" A = convert(Matrix{T}, triu(randn(n,n)))\n",
" x = A \\ b\n",
" c[i] = max(log10(abs(norm(A*x)-norm(b))/\u03b5), -2) #arbitrary floor for plotting 0 on log scale\n",
" resid = norm(abs(b - A*x))\n",
" rhs = \u03b3 * norm(A) * norm(x)\n",
" r[i] = resid/rhs\n",
"end"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 13
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"plot(x=c, y=r, Guide.XLabel(\"log10(c)\"), Guide.YLabel(\"lhs/rhs\"))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"\n",
"\n"
],
"metadata": {},
"output_type": "display_data",
"text": [
"D3(120.0,80.0,IOBuffer([0x66,0x75,0x6e,0x63,0x74,0x69,0x6f,0x6e,0x20,0x64 \u2026 0x20,0x3d,0x20,0x64,0x72,0x61,0x77,0x0a,0x7d,0x0a],true,true,true,false,382058,9223372036854775807,382059),0,[],[],1,3,[0xebf397648d57c064=>([56.75578406500111,56.354913590113426,52.69882990819332,54.183497079931854,54.005450874262856,51.868729398992976,54.253616289692864,52.33910055295948,57.6352758451944,50.14191392606996 \u2026 53.335518793023176,56.36651635781103,54.95243050136375,56.99367490040932,47.54277499332237,57.493080370459,52.79221485928805,53.36356047597386,56.935877071024585,48.893016270601166],1),0x7c2384964ac76785=>([53.80661954144316,29.789849056603774,51.3162495056162,64.28971607041647,51.3162495056162,59.58907969429545,61.89516089946106,29.789849056603774,69.5505598376033,58.787359613097095 \u2026 53.80661954144316,51.3162495056162,62.58159265618741,63.87865746936825,65.1492311431727,51.3162495056162,56.296989577270125,57.09870965846848,57.09870965846848,51.3162495056162],0)],true,false,nothing,true)"
]
},
{
"html": [],
"metadata": {},
"output_type": "pyout",
"prompt_number": 14,
"text": [
"Plot(...)"
]
}
],
"prompt_number": 14
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The original Julia test doesn't seem to correlate with the test derived from backward error result either. But at least none of the tests fail so far."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"x, y = hist(r, 25)\n",
"plot(x=x, y=y/(t*(x[2]-x[1])), Geom.bar, Guide.XLabel(\"lhs/rhs\"), Guide.YLabel(\"Density\"))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"\n",
"\n"
],
"metadata": {},
"output_type": "display_data",
"text": [
"D3(120.0,80.0,IOBuffer([0x66,0x75,0x6e,0x63,0x74,0x69,0x6f,0x6e,0x20,0x64 \u2026 0x20,0x3d,0x20,0x64,0x72,0x61,0x77,0x0a,0x7d,0x0a],true,true,true,false,20564,9223372036854775807,20565),0,[],[],0,3,Dict{Uint64,(Any,Int64)}(),true,false,nothing,true)"
]
},
{
"html": [],
"metadata": {},
"output_type": "pyout",
"prompt_number": 15,
"text": [
"Plot(...)"
]
}
],
"prompt_number": 15
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"mean(r), std(r)/sqrt(t-1)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 16,
"text": [
"(0.005644636805875966,4.95668601105437e-5)"
]
}
],
"prompt_number": 16
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Summary\n",
"\n",
"This Issue Note describes two replacement tests for the existing test for the linear solve with a triangular matrix. Numerical testing indicates that the new tests are more robust and do not exhibit any spurious failures when applied to a moderate number of randomly sampled matrices."
]
}
],
"metadata": {}
}
]
}