title: Julia: isposdef() fails for large matrices?
body: <p>I have a positive definite covariance matrix <code>C</code> of size 3n x 3n constructed from n^2  blocks of size 3x3.</p>
  <p>Running <code>MvNormal</code> (with e.g a zero mean vector) on this matrix to draw Gaussian random vectors, I am getting the error</p>
  <pre><code>PosDefException: matrix is not positive definite; Cholesky factorization failed.
  </code></pre>
  <p>and indeed checking <code>isposef(C)</code> returns <code>false</code> when n becomes too large. However my matrix should be positive definite for any n, so it seems that there is some kind of numerical instability (perhaps due to the determinant becoming too small or too large beyond machine precision).</p>
  <p>The reproducible code I am using to generate C is below:</p>
  <pre><code>#######################################################
  # inputs
  grid_size=10
  l_sq = 1 
  xmax = 2
  #######################################################
  
  # kernel function used to construct covariance matrix C
  function corr(x, y, l_sq)
      v=x-y
      d_sq=sum(v.^2)
  
      n = d_sq/(2l_sq)
      return exp(-n)*(Matrix{Float64}(I, length(x), length(x)))
  end 
  
  nb_grid_points = grid_size^3
  gaussian_vector_dim = 3*nb_grid_points
  
  oneD_grid = LinRange(-xmax, xmax, grid_size)
  
  # get input set X which indexes grid points
  threeD_grid = collect.(Iterators.product(oneD_grid, oneD_grid, oneD_grid))
  grid_points = vec(reshape(threeD_grid,:,1))
  
  ########################################
  # build C by blocks
  C = Array{Float64}(undef, gaussian_vector_dim, gaussian_vector_dim)
  
  for i in 1:nb_grid_points
      for j in 1:nb_grid_points
          #block covariance matrix C consist of DxD correlation-function matrices K_i,j for i,j=1,...,nb_grid_points
          C[3*(i-1)+1:(3*i),(3*(j-1)+1):(3*j)] = corr(grid_points[i], grid_points[j],l_sq)
      end
  end
  
  #########################################
  # plot covariance matrix
  plt.imshow(C,cmap=&quot;Blues&quot;, interpolation=&quot;none&quot;) 
  plt.colorbar() 
  plt.title(&quot;Covariance matrix&quot;)
  #########################################
  print(&quot;C is symmetric:&quot;,issymmetric(C))
  print(&quot;\ndet C=&quot;,det(C))
  print(&quot;\nC is positive definite=&quot;,isposdef(C))
  </code></pre>
  <p>Maintaining,     <code>l_sq = 1 , xmax = 2</code>, the code above gives <code>isposdef(C) = false</code> when    <code>grid_size=10</code> but  <code>isposdef(C)=true</code> if <code>grid_size</code> is 9 or less.</p>
  <p>Why is this failure occurring and how can I fix it? Perhaps I can help Julia by indicating that the covariance matrix is sparse?</p>