Solve vs. LinearSolve - ill solutions using both methods

Question

I'm wishing you a nice day. My question is related to this question previously asked by me, answered by PlatoManiac. Although I found a bug in his post I decided to leave his answer accepted and ask about the issue separately, because the bug is somewhat unrelated to what I'd trying to achieve in the aforementioned question (and it was not his job to dig through the fits he calculated to see the problem).

The problem is with functions Solve and LinearSolve, I'll try to explain. Please follow these commands exactly:

Needs["NDSolve`FEM`"]
region = ImplicitRegion[
   x^2 + y^2 - (1/2)^2 >= 0, {{x, -1, 1}, {y, -1, 1}}];
mesh = ToElementMesh[region, "MaxBoundaryCellMeasure" -> 0.05, 
   "MaxCellMeasure" -> 0.005, "MeshOrder" -> 1];
grid = mesh[[1]];
len = Length[grid]
765

For every point I'll find 6 nearest points:

neighbour = Table[Nearest[grid, grid[[i]], 6], {i, 1, len}]

So for example the first point has following neighbours:

neighbour[[1]]
{{0.5, -1.11022*10^-15}, {0.497508, -0.04986}, {0.497508, 0.04986}, {0.54888, 0.0274355}, {0.557905, -0.0278866}, {0.490056, -0.0992229}}

Now let's try the very same fit as in the answer:

fit[x_, y_] := a1 + a2 x + a3 y + a4 x y + a5 x^2 + a6 y^2

Command fit@@@neighbour[[i]] will map function fit to the six neighbour points of the i-th point of the grid. I want to solve equation of the form:

Solve[fit@@@neighbour[[i]] == {a,b,c,d,e,f}, {a1,a2,a3,a4,a5,a6}]

for some i. There are several equations which has apparently no solution:

Position[Table[
  Solve[fit@@@neighbour[[i]] == {a, b, c, d, e, f}, {a1, a2, a3, a4,
     a5, a6}], {i, 1, len}], {}]
{{1}, {33}, {85}, {126}}

So, let's look at them closely; first, make an matrix of every one of them:

matrix1 = 
 Normal@CoefficientArrays[
    fit@@@neighbour[[1]], {a1, a2, a3, a4, a5, a6}][[2]]; matrix1 // MatrixForm

$$ \tiny \left( \begin{array}{cccccc} 1. & 0.5 & -\text{1.1102230246067903$\grave{ }$*${}^{\wedge}$-15} & -\text{5.551115123033953$\grave{ }$*${}^{\wedge}$-16} & 0.25 & \text{1.2325951643670497$\grave{ }$*${}^{\wedge}$-30} \\ 1. & 0.497508 & -0.04986 & -0.0248057 & 0.247514 & 0.00248602 \\ 1. & 0.497508 & 0.04986 & 0.0248057 & 0.247514 & 0.00248602 \\ 1. & 0.54888 & 0.0274355 & 0.0150588 & 0.301269 & 0.000752708 \\ 1. & 0.557905 & -0.0278866 & -0.0155581 & 0.311258 & 0.000777663 \\ 1. & 0.490056 & -0.0992229 & -0.0486248 & 0.240155 & 0.00984519 \\ \end{array} \right)$$

matrix33 = 
     Normal@CoefficientArrays[
        fit@@@neighbour[[33]], {a1, a2, a3, a4, a5, a6}][[2]];
matrix85 = 
     Normal@CoefficientArrays[
        fit@@@neighbour[[85]], {a1, a2, a3, a4, a5, a6}][[2]];
matrix126 = 
     Normal@CoefficientArrays[
        fit@@@neighbour[[126]], {a1, a2, a3, a4, a5, a6}][[2]];

And let's try LinearSolve instead of Solve:

LinearSolve[matrix1, {a, b, c, d, e, f}]

And the solution indeed exists!

{4247.93 (1. a - 1.03612 b - 0.282395 c - 0.0880603 d + 0.0728978 e + 
0.333914 f), -16167.9 (1. a - 1.03725 b - 0.282848 c - 0.0875778 d + 
0.0734418 e + 0.334234 f), -1183.67 (1. a - 0.974568 b - 0.375339 c + 
0.0597934 d - 0.0501421 e + 0.340256 f), 2379.2 (1. a - 0.98304 b - 
0.366867 c + 0.0597934 d - 0.0501421 e + 0.340256 f), 15348.1 (1. a - 
1.03823 b - 0.283277 c - 0.0870207 d + 0.074025 e + 0.3345 f), 
-1262.48 (1. a - 0.854177 b - 0.34683 c - 0.0664458 d + 0.0429523 e + 0.2245 f)}

The same it's with other three cases when Solve found no solution - LinearSolve found solution for every grid point, which is good:

sols = Table[
      LinearSolve[
       Normal@CoefficientArrays[
          fit@@@neighbour[[1]], {a1, a2, a3, a4, a5, a6}][[2]], {a, b, 
        c, d, e, f}], {i, 1, len}];
Position[sols, {}]
{}

But I sense another problem now! Let's find some enormously big coefficients in sols:

Position[sols, _?(# > 10^10 &)]

The first such solution is sols[[79]] and MON DIEU look at that!

sols[[79]]
{2.08455*10^44 a - 2.08455*10^44 b - 6.94849*10^43 c - 2046.98 d + 
  1678.53 e + 6.94849*10^43 f, 
 4.43182*10^44 a - 4.43182*10^44 b - 1.47727*10^44 c - 4212.86 d + 
  3435.28 e + 1.47727*10^44 f, -4.86799*10^43 a + 4.86799*10^43 b + 
  1.62266*10^43 c - 4.9869*10^-14 d + 4.9869*10^-14 e - 
  1.62266*10^43 f, -4.86799*10^43 a + 4.86799*10^43 b + 
  1.62266*10^43 c - 4.9869*10^-14 d + 4.9869*10^-14 e - 
  1.62266*10^43 f, 
 2.34727*10^44 a - 2.34727*10^44 b - 7.82423*10^43 c - 2165.88 d + 
  1756.75 e + 7.82423*10^43 f, -4.27219*10^16 a + 4.27219*10^16 b + 
  1.42406*10^16 c - 1.42406*10^16 f}

Would you look at that...10 to the 44...that is REALLY big number! Let's try Solve again and solve 79-th equation:

Solve[fit @@@ neighbour[[79]] == {a, b, c, d, e, f}, {a1, a2, a3, a4, 
  a5, a6}]
{{a1 -> 1.81587*10^16 a - 1.81587*10^16 b - 6.05291*10^15 c - 
    292.579 d + 260.359 e + 6.05291*10^15 f, 
  a2 -> 3.8606*10^16 a - 3.8606*10^16 b - 1.28687*10^16 c - 
    482.94 d + 420.2 e + 1.28687*10^16 f, 
  a3 -> -4.24056*10^15 a + 4.24056*10^15 b + 1.41352*10^15 c - 
    409.701 d + 331.182 e - 1.41352*10^15 f, 
  a4 -> -4.24056*10^15 a + 4.24056*10^15 b + 1.41352*10^15 c - 
    409.701 d + 331.182 e - 1.41352*10^15 f, 
  a5 -> 2.04473*10^16 a - 2.04473*10^16 b - 6.81578*10^15 c - 
    190.361 d + 159.842 e + 6.81578*10^15 f, 
  a6 -> 440.433 a - 640.433 b - 80.1445 c - 8.79132*10^-12 d + 
    7.51383*10^-12 e + 280.144 f}}

Clearly, this solution is absolutely different from the one obtained by LinearSolve. So now there are several problems: at least one of the solutions obtained by Solve and LinearSolve must be wrong, both of them can't be correct. Then I can't use Solve to create set of interpolants at each grid point, because at some points, no solution exists ({}), but I can't use LinearSolve neither, because I'm not sure if such monstrous numbers (10^44) as a result are correct...I need this to approximate a derivatives in a grid points (a,b,c,d,e,f stands for some yet unknown functional values) and 30 orders of magnitude is really huge error...so any help would be appreciated. Thanks in advance!

P.S.: each set of the six nearest grid points to the every grid point seem to be well distributed to me...neither two points are too close nor too far away. Try ListPlot[neighbour[[i]]] to check that. I'm assuming that the solutions to the discussed equations must have somewhat "normal" numerical coefficients in them, not the mad things like 10^16 or 10^44...

I think that is why FindRoot does not converge in this other question of mine - let's imagine the Newton's method: algorithm calculates tangent line, makes a movement but the value 10^44 (even 10^16) is SO HUGE that in the next step the position is too far away from solution.