fix #535 for GaussNewton algorithm failure for sparse jacobian matrices

[Iskaban10]

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

This PR addresses Issue #535, where GaussNewton (via NewtonDescent) diverges numerically when provided with a sparse jac_prototype in a NonlinearLeastSquaresProblem.
Add any other context about the problem here.
The root cause was that newtondescent skipped the Normal Form ($J^T J \Delta u = -J^T f$) for sparse matrices if the linear solver didn't explicitly signal a need for square matrices.

[Iskaban10]

Is this approach correct? I have used SparseArray.issparse() to check for sparse Jacobians to force them to be squared for normal forms.

[oscardssmith]

why do we want to use the normal form for sparse matrices?

[ChrisRackauckas]

normal form should be less numerically stable, just slightly faster

[Iskaban10]

why do we want to use the normal form for sparse matrices?
Yeah I thought the error was occuring because the sparse matrix was in a rectangular form making the system ill-defined. Therefore attempted a brute force technique but it isn't the right solution.

[Iskaban10]

normal form should be less numerically stable, just slightly faster

Based on my understnding of the issue till now, I think the best way here is to apply a line search technique so that the step size is shrinked in case such high residual norms appear.
Would highly appreciate any advice.

[ChrisRackauckas]

Then choose a method with a line search, or LM?

[Iskaban10]

Then choose a method with a line search, or LM?

Thanks for the guidance. After digging a bit deeper (along with help from AI), I found out that in newton.jl, in case of sparse matrices, the in-place multiplication was silently dropping entries that fell outside the sparsity pattern of the initialized cache, as u changed. This resulted in a corrupted normal equations matrix and the subsequent garbage values.
I have made a new PR based on this.