This is a linear model. Therefore the Hessians computed later should agree with the Gauss-Newton approximation.

xs=Table[i/10, {i,0,9}];
noise=10^-5;
yos=MapIndexed[Function[{x,i},2x^2(1-x)+noise (-1)^(First[i]+1)],xs];
yc[a_,b_,c_,x_]:=-x^3+a x^2+b x + c
residuals[a_,b_,c_,K_]:=MapThread[Function[{x,yo}, yo-K yc[a,b,c,x]],{xs,yos}]
L[a_,b_,c_,K_]:=Simplify[1/2Plus@@(residuals[a,b,c,K]^2)]

vals={a->0.5,b->0.3,c->0.2}

J[a,b,c]:=Table[Map[Function[x,D[yc[a,b,c,x],α]],xs],{α,{a,b,c}}]
J[a,b,c].J[a,b,c]\[Transpose]/.vals//N//MatrixForm

diffKstar[a,b,c]:=Table[D[Kstar[a,b,c],α],{α,{a,b,c}}]
diffKstar[a,b,c]/.vals//N

{Lstar[a_,b_,c_],ruleKstar[a_,b_,c_]}=Minimize[L[a,b,c,K],K]//Simplify;
Kstar[a_,b_,c_]=K/.ruleKstar[a,b,c];
Hgn[a_,b_,c_]:=Table[D[residuals[a,b,c,Kstar[a,b,c]],α].D[residuals[a,b,c,Kstar[a,b,c]],β],{α,{a,b,c}},{β,{a,b,c}}]
bgn[a_,b_,c_]:=Table[-residuals[a,b,c,Kstar[a,b,c]].D[residuals[a,b,c,Kstar[a,b,c]],α],{α,{a,b,c}}]

Kstar[a,b,c]/.vals//InputForm

L[a,b,c,Kstar[a,b,c]]/.vals//InputForm

AA=Hgn[a,b,c]/.vals;
AA//MatrixForm
AA//InputForm

bb=bgn[a,b,c]/.vals;
bb//InputForm

{λ,x}=Eigensystem[AA];
λ
x//MatrixForm

U=CholeskyDecomposition[AA];
U//InputForm
U//MatrixForm

s=LinearSolve[AA,bb];
s//InputForm

Now let's add a quadratic penalty...

q[a_,b_,c_]:=1/2(a-b+c-2)^2
L1[a_,b_,c_]:=Lstar[a,b,c]+q[a,b,c]
L1[a,b,c]/.vals//InputForm

b1gn[a_,b_,c_]:=bgn[a,b,c]-D[q[a,b,c],{{a,b,c}}]
b1gn[a,b,c]/.vals//InputForm

H1gn[a_,b_,c_]:=Hgn[a,b,c]+D[q[a,b,c],{{a,b,c}},{{a,b,c}}]
H1gn[a,b,c]/.vals//InputForm

This is a linear model. Therefore the Hessians computed later should agree with the Gauss - Newton approximation.

ClearAll[x,a,L]
noise=10^-5
ts=Table[i/10, {i,0,9}];
f[t_]:=t-t^2+2t^3
yos=MapIndexed[Function[{t,i},f[t]+f[1-t]+f[-t]+noise(-1)^(First[i]+1)],ts];
yc[x_,t_]:=Sum[x[[i]]t^i,{i,4}]
residuals[x_,a_]:=MapThread[Function[{t,yo},(yo-a[[1]] yc[x,t]-a[[2]] yc[x,1-t]-a[[3]] yc[x,-t])^2],{ts,yos}]
L[x_,a_]:=Plus@@residuals[x,a]
xx={x1,x2,x3,x4}
aa={a1,a2,a3}
L[xx,aa]//Simplify;

vals={x1->1,x2->2,x3->3,x4->4}

ClearAll[AA,bb]
Yc[x_]:=Map[Function[t,{yc[x,t],yc[x,1-t],yc[x,-t]}],ts]
AA=Yc[xx]\[Transpose].Yc[xx];
bb=Yc[xx]\[Transpose].yos;
N[AA/.vals,15]//MatrixForm
N[bb/.vals,15]//MatrixForm
astar=LinearSolve[AA,bb]//Simplify;
U=CholeskyDecomposition[AA];
N[astar/.vals,15]
astarrule={a1->astar[[1]],a2->astar[[2]],a3->astar[[3]]};
{D[L[xx,aa],a1],D[L[xx,aa],a2],D[L[xx,aa],a3]}/.astarrule/.vals

Hs=N[Table[D[L[xx,aa]/.astarrule,u,v],{u,{x1,x2,x3,x4}},{v,{x1,x2,x3,x4}}]/.vals,15]/2
Hs//MatrixForm

f=-N[Table[D[L[xx,aa]/.astarrule,u],{u,{x1,x2,x3,x4}}]/.vals,15]/2

index[i_,j_,n_]:=(i-1)n-(i-1)i/2+j
PackedUAsMat[u_,i_,j_,n_]:=If[i>j,u[[index[j,i,n]]],u[[index[i,j,n]]]]
SymmetricMatrixFromUpperDiagonal[u_]:=Block[{n},
n=(Sqrt[1+8Length[u]]-1)/2;
Table[PackedUAsMat[u,i,j,n],{i,1n},{j,1,n}]]