We currently solve, J^TJ = -J^T f in parallel with partial pivoting of Gaussian elimination. Since we know that the left hand side is positive semi definite, one does not need to do Gaussian elimination with partial pivoting. We could instead use Cholesky decomposition.