<p dir="auto">I wrote a function to rotate the coordinates of a 2D MvNormal distribution, and in the course of testing it I discovered the Cholesky factorization throws an error depending on the input with bit-level sensitivity. Is this a bug? For my use case, the workaround of wrapping the input in <code class="notranslate">Hermitian()</code> doesn't work, because I can't pass a <code class="notranslate">Hermitian</code> to the <code class="notranslate">MvNormal</code> constructor, which is calling <code class="notranslate">cholesky</code> for me.</p>
<p dir="auto">Example:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="julia&gt; VERSION
v&quot;1.1.0&quot;

julia&gt; rot1 = [1 -1; 1 1] / sqrt(2)
2×2 Array{Float64,2}:
 0.707107  -0.707107
 0.707107   0.707107

julia&gt; rot2 = [cos(π/4) -sin(π/4); sin(π/4) cos(π/4)]
2×2 Array{Float64,2}:
 0.707107  -0.707107
 0.707107   0.707107

julia&gt; cholesky(rot1 * [25.0 0.0; 0.0 4.0] * rot1')
Cholesky{Float64,Array{Float64,2}}
U factor:
2×2 UpperTriangular{Float64,Array{Float64,2}}:
 3.80789  2.75744
  ⋅       2.62613

julia&gt; cholesky(rot2 * [25.0 0.0; 0.0 4.0] * rot2')
ERROR: PosDefException: matrix is not Hermitian; Cholesky factorization failed.
Stacktrace:
 [1] checkpositivedefinite(::Int64) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/factorization.jl:11
 [2] #cholesky!#97(::Bool, ::Function, ::Array{Float64,2}, ::Val{false}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/cholesky.jl:182
 [3] #cholesky#101 at ./none:0 [inlined]
 [4] cholesky at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/cholesky.jl:275 [inlined] (repeats 2 times)
 [5] top-level scope at none:0"><pre class="notranslate"><code class="notranslate">julia&gt; VERSION
v"1.1.0"

julia&gt; rot1 = [1 -1; 1 1] / sqrt(2)
2×2 Array{Float64,2}:
 0.707107  -0.707107
 0.707107   0.707107

julia&gt; rot2 = [cos(π/4) -sin(π/4); sin(π/4) cos(π/4)]
2×2 Array{Float64,2}:
 0.707107  -0.707107
 0.707107   0.707107

julia&gt; cholesky(rot1 * [25.0 0.0; 0.0 4.0] * rot1')
Cholesky{Float64,Array{Float64,2}}
U factor:
2×2 UpperTriangular{Float64,Array{Float64,2}}:
 3.80789  2.75744
  ⋅       2.62613

julia&gt; cholesky(rot2 * [25.0 0.0; 0.0 4.0] * rot2')
ERROR: PosDefException: matrix is not Hermitian; Cholesky factorization failed.
Stacktrace:
 [1] checkpositivedefinite(::Int64) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/factorization.jl:11
 [2] #cholesky!#97(::Bool, ::Function, ::Array{Float64,2}, ::Val{false}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/cholesky.jl:182
 [3] #cholesky#101 at ./none:0 [inlined]
 [4] cholesky at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/cholesky.jl:275 [inlined] (repeats 2 times)
 [5] top-level scope at none:0
</code></pre></div>
<p dir="auto">Related issues from Discourse:</p>
<ul dir="auto">
<li><a href="https://discourse.julialang.org/t/multivariate-normal-with-positive-semi-definite-covariance-matrix/3029" rel="nofollow">https://discourse.julialang.org/t/multivariate-normal-with-positive-semi-definite-covariance-matrix/3029</a></li>
<li><a href="https://discourse.julialang.org/t/is-this-a-bug-with-cholesky/16970" rel="nofollow">https://discourse.julialang.org/t/is-this-a-bug-with-cholesky/16970</a></li>
<li><a href="https://discourse.julialang.org/t/numerical-stability-of-cholesky-factorization/20063" rel="nofollow">https://discourse.julialang.org/t/numerical-stability-of-cholesky-factorization/20063</a></li>
</ul>
<p dir="auto">I have a function where i determine the hessian with finite difference using the <code class="notranslate">hessian</code> function in the Calculus package. Then I test if the hessian is positive definite with <code class="notranslate">isposdef</code> this always return false because, <code class="notranslate">ishermitian</code> is called at: </p><div class="Box Box--condensed my-2">
  <div class="Box-header f6">
    <p class="mb-0 text-bold">
      <a href="https://github.com/JuliaLang/julia/blob/5e2a6173eca65e78f787a3525bae18857d283229/base/linalg/dense.jl#L27-L28">julia/base/linalg/dense.jl</a>
    </p>
    <p class="mb-0 color-fg-muted">
        Lines 27 to 28
      in
      <a data-pjax="true" class="commit-tease-sha" href="/JuliaLang/julia/commit/5e2a6173eca65e78f787a3525bae18857d283229">5e2a617</a>
    </p>
  </div>
  <div itemprop="text" class="Box-body p-0 blob-wrapper blob-wrapper-embedded data">
    <table class="highlight tab-size mb-0 js-file-line-container" data-tab-size="8" data-paste-markdown-skip="">

        <tbody><tr class="border-0">
          <td id="L27" class="blob-num border-0 px-3 py-0 color-bg-default" data-line-number="27"></td>
          <td id="LC27" class="blob-code border-0 px-3 py-0 color-bg-default blob-code-inner js-file-line"> <span class="pl-en">isposdef!</span><span class="pl-c1">{T&lt;:BlasFloat}</span>(A<span class="pl-k">::</span><span class="pl-c1">StridedMatrix{T}</span>, UL<span class="pl-k">::</span><span class="pl-c1">Symbol</span>) <span class="pl-k">=</span> LAPACK<span class="pl-k">.</span><span class="pl-c1">potrf!</span>(<span class="pl-c1">char_uplo</span>(UL), A)[<span class="pl-c1">2</span>] <span class="pl-k">==</span> <span class="pl-c1">0</span> </td>
        </tr>

        <tr class="border-0">
          <td id="L28" class="blob-num border-0 px-3 py-0 color-bg-default" data-line-number="28"></td>
          <td id="LC28" class="blob-code border-0 px-3 py-0 color-bg-default blob-code-inner js-file-line"> <span class="pl-en">isposdef!</span>(A<span class="pl-k">::</span><span class="pl-c1">StridedMatrix</span>) <span class="pl-k">=</span> <span class="pl-c1">ishermitian</span>(A) <span class="pl-k">&amp;&amp;</span> <span class="pl-c1">isposdef!</span>(A, <span class="pl-c1">:U</span>) </td>
        </tr>
    </tbody></table>
  </div>
</div>
<br>
before the <code class="notranslate">LAPACK.potrf!</code> is called. Due to numerical inaccuracy will <code class="notranslate">ishermitian</code> always return false because <code class="notranslate">A[i,j] != ctranspose(A[j,i])</code> in <div class="Box Box--condensed my-2">
  <div class="Box-header f6">
    <p class="mb-0 text-bold">
      <a href="https://github.com/JuliaLang/julia/blob/5e2a6173eca65e78f787a3525bae18857d283229/base/linalg/generic.jl#L282">julia/base/linalg/generic.jl</a>
    </p>
    <p class="mb-0 color-fg-muted">
         Line 282
      in
      <a data-pjax="true" class="commit-tease-sha" href="/JuliaLang/julia/commit/5e2a6173eca65e78f787a3525bae18857d283229">5e2a617</a>
    </p>
  </div>
  <div itemprop="text" class="Box-body p-0 blob-wrapper blob-wrapper-embedded data">
    <table class="highlight tab-size mb-0 js-file-line-container" data-tab-size="8" data-paste-markdown-skip="">

        <tbody><tr class="border-0">
          <td id="L282" class="blob-num border-0 px-3 py-0 color-bg-default" data-line-number="282"></td>
          <td id="LC282" class="blob-code border-0 px-3 py-0 color-bg-default blob-code-inner js-file-line"> <span class="pl-k">if</span> A[i,j] <span class="pl-k">!=</span> <span class="pl-c1">ctranspose</span>(A[j,i]) </td>
        </tr>
    </tbody></table>
  </div>
</div>
, will return true for some <code class="notranslate">i,j</code> due to numerical inaccuracy. If I change the line to <code class="notranslate">!isapprox(A[i,j],ctranspose(Aj,i]))</code> the <code class="notranslate">ishermitian</code> returns true and hence also call the <code class="notranslate">LAPACK.potrf!</code> to check if it is positive definite.<p></p>
<p dir="auto">So from my point of view, should <code class="notranslate">ishermitian</code> test approximately equality.</p>