МНК с особой точкой

[pvyu]

Код для обращения таков (взят отсюда):

// Обращение квадратной матрицы
template<typename TT>
inline static bool InvertMN(int N, TT **M,
                            int         *GausOstatok=NULL,
                            long double *GausDet =   NULL,
                            long double *GausMinved= NULL) {
  // (c) Drobotenko http://www.drobotenko.com/bib_rus.html

  int gaus_ostatok;               // дефект линейного преобразования  0 == ОК
  long double gaus_deter;     // определитель в случае успешного обращения
  long double gaus_minved;  // минимальный ведущий элемент
                                           // можно использовать для оценки точности
  int *rn, *cn;
  int j, k;

  #define _NaN() (double&)*"Implementation dependent";
  //  заполнитель для неразрешенных измерений
  //  при невозможности обращения

  rn = new int [N];
  cn = new int [N];
  for(j=N; j--; )
    rn[j]=cn[j]=j;
  gaus_minved=1e99;
  gaus_deter=1;
  for(gaus_ostatok=N; gaus_ostatok; gaus_ostatok--) {
    int jved,kved;
    long double vved=-1, t;

    for(j=N; j--; ) {  // поиск ведущего
       if(~rn[j])
         for(k=N; k--; )
           if(~cn[k])
             if(vved<fabsl(M[j][k]))
                vved=fabsl(M[j][k]),jved=j,kved=k;  }

    if(gaus_minved>vved)
      gaus_minved=vved;
    gaus_deter*=M[jved][kved];
    if(vved<1e-99) {          // Для long double можно уменьшить
      for(j=N; j--; ) {
         if(~rn[j])
           for(k=N; k--; )
             M[j][k]=_NaN();  // ?
         if(~cn[j])
           for(k=N; k--; )
             M[k][j]=_NaN();  }
      delete[] rn;
      delete[] cn;
      return false;  }
    int jt=rn[jved],kt=cn[kved];
    // перестановки
    for(j=N; j--; )
      t=M[kt][j],M[kt][j]=M[jved][j],M[jved][j]=t;
    for(j=N; j--; )
      t=M[j][jt],M[j][jt]=M[j][kved],M[j][kved]=t;

    rn[jved]=rn[kt];
    cn[kved]=cn[jt];
    rn[kt]=cn[jt]=-1;

    vved=M[kt][jt];   M[kt][jt]=1;
    for(j=N; j--; )  {
      if(j==kt)
        continue;
      long double mul=M[j][jt]/vved;
      M[j][jt]=0;
      for(k=N; k--; )
        M[j][k]-=M[kt][k]*mul;  }  // самый внутренний цикл ровно N^3 операций
    for(k=N; k--; )
      M[kt][k]/=vved;
  }
  if (GausOstatok) *GausOstatok=gaus_ostatok;
  if (GausDet)        *GausDet =   gaus_deter;
  if (GausMinved)   *GausMinved= gaus_minved;

  delete[] rn;
  delete[] cn;
  return true;
}

Код работоспособный, при перемножении матриц $A^{-1}\cdot A$ размером 3x3 погрешности имеют порядок -19.
Матрицы, используемые в выражениях, заданы, предположительно, верно.
Выражение $(X^TX)^{-1} $ в случае $X=\begin{bmatrix}5 & 20 \\20 & 90 \end{bmatrix}$ вычисляется так:
$(X^TX)^{-1}=\begin{bmatrix}3,4 & -0,76 \\-0,76 & 0,17 \end{bmatrix}$,
аналогичные вычисления в MatLab дают те же результаты.