Решил выложить пару файлов с реализацией линейного метода наименьших квадратов. Может быть, кому-нибудь пригодится.
Буду рад вопросам, замечаниям и предложениям...

---

Небольшой класс для представления матриц:

// file Matrix.h

#ifndef MATRIX_H_INCLUDED
#define MATRIX_H_INCLUDED

#include <vector>
#include <istream>
#include <ostream>
#include <string>
#include <sstream>
#include <iterator>
#include <iomanip>

template<class T>
class Matrix
{
public:
    Matrix() {}
    Matrix(int rows, int cols) : rows_(rows), cols_(cols), data(rows_ * cols_) {}

    int Rows() const { return rows_; }
    int Cols() const { return cols_; }

    void Resize(int rows, int cols)
    {
        rows_ = rows; cols_ = cols;
        data.resize(rows_ * cols_);
    }
    
    T& operator()(int row, int col) { return data[row * cols_ + col]; }
    T operator()(int row, int col) const { return data[row * cols_ + col]; }
    
    T* Row(int row) { return &data[row * cols_]; }
    const T* Row(int row) const { return &data[row * cols_]; }

    void Read(std::istream& is)
    {
        data.clear();
        rows_ = 0; cols_ = 0;
        int old_size;
        
        std::string s;
        while (std::getline(is, s))
        {
            std::istringstream iss(s.c_str());
            std::copy(std::istream_iterator<double>(iss),
                std::istream_iterator<double>(),
                std::back_inserter(data));
            
            int n = data.size();
            if (cols_ == 0)
                cols_ = n;
            else
                if ((n - old_size) != cols_)
                {
                    data.clear();
                    rows_ = 0;
                    cols_ = 0;
                    throw std::logic_error("invalid matrix");
                }
                
            old_size = n;
            ++rows_;
        }
    }

    void Write(std::ostream& os) const
    {
        os.flags(std::ios::fixed);
        for (int i=0; i<rows_; ++i)
        {
            for (int j=0; j<cols_; ++j)
                os << std::setw(15) << std::setprecision(8) << (*this)(i, j);
            os << '\n';
        }
    }

    Matrix Inverse() const
    {
        if (rows_ != cols_)
            throw std::invalid_argument("trying to invert non-square matrix");
        
        Matrix A(rows_, cols_);
        std::fill(A.data.begin(), A.data.end(), 0);
        for (int i=0; i<rows_; ++i)
            A(i,i) = 1;
        
        Matrix B(*this);

        for (i=0; i<rows_; ++i)
        {
            int j;
            for (j=i; (j<rows_) && (B(i,j) == 0); ++j);
            
            if (j == rows_)
                throw std::logic_error("invalid matrix");
            
            if (i != j)
                for (int k=0; k<rows_; ++k)
                {
                    std::swap(A(i,k), A(j,k));
                    std::swap(B(i,k), B(j,k));
                }
                
                double d = B(i,i);
                for (j=0; j<rows_; ++j)
                {
                    A(i,j) /= d;
                    B(i,j) /= d;
                }
                
                for (j=i+1; j<rows_; ++j)
                {
                    d = B(j,i);
                    for (int k=0; k<rows_; ++k)
                    {
                        A(j,k) -= A(i,k) * d;
                        B(j,k) -= B(i,k) * d;
                    }
                }
        }
        
        for (i=rows_-1; i>0; --i)
            for (int j=0; j<i; ++j)
            {
                double d = B(j,i);
                for (int k=0; k<rows_; ++k)
                {
                    A(j,k) -= A(i,k) * d;
                    B(j,k) -= B(i,k) * d;
                }
            }
            
        return A;
    }
    
    friend Matrix operator*(const Matrix&, const Matrix&);

protected:
    int rows_, cols_;
    std::vector<T> data;
};

template<class T>
Matrix<T> operator*(const Matrix<T>& A, const Matrix<T>& B)
{
    if (A.Cols() != B.Rows())
        throw std::invalid_argument("invalid matrix sizes");

    Matrix<T> C;

    C.Resize(A.Rows(), B.Cols());

    std::fill(C.data.begin(), C.data.end(), 0);

    for (int i=0; i<C.Rows(); ++i)
        for (int j=0; j<C.Cols(); ++j)
            for (int k=0; k<A.Cols(); ++k)
                C(i, j) += A(i, k) * B(k, j);
    
    return C;
}

#endif // MATRIX_H_INCLUDED

Собственно, сам МНК:

// file LinFit.h

#ifndef LINFIT_H_INCLUDED
#define LINFIT_H_INCLUDED

#include "Matrix.h"

template<class T, class Model>
Matrix<double> LinearFit(const Matrix<T> &D, const Model &f)
{
    Matrix<double> B(f.Dim(), 1);
    Matrix<double> P(f.Dim(), f.Dim());

    for (int i=0; i<f.Dim(); ++i)
    {
        B(i, 0) = 0;
        for (int j=0; j<D.Rows(); ++j)
            B(i, 0) += D(j, 2) * f(i, D.Row(j));
            
        for (j=0; j<f.Dim(); ++j)
        {
            P(i, j) = 0;
            for (int k=0; k<D.Rows(); ++k)
                P(i, j) += f(i, D.Row(k)) * f(j, D.Row(k));
        }
    }
    
    P = P.Inverse();
    
    return P * B;
}

#endif // LINFIT_H_INCLUDED

Пример использования:

// file PolyFit.cpp

#include "matrix.h"
#include "LinFit.h"

#include <iostream>
#include <fstream>

using std::cout;
using std::ios;

struct Model
{
    int operator()(int k, const int *x) const
    {
        switch(k)
        {
        case 0: return 1;
        case 1: return x[0];
        case 2: return x[1];
        case 3: return x[0]*x[0];
        case 4: return x[1]*x[1];
        case 5: return x[0]*x[1];
        default: return 0;
        }
    }

    int Dim() const
    {
        return 6;
    }
};

int main(int argc, char *argv[])
{
    if (argc < 2)
    {
        cout << "USAGE: polyfit <infile>\n";
        return 0;
    }

    Matrix<int> D;
    std::ifstream fin(argv[1]);
    D.Read(fin);
    
    Matrix<double> B;
    
    try
    {
        B = LinearFit(D, Model());
    }
    catch(std::exception &e)
    {
        cout << e.what() << '\n';
    }
    
    B.Write(cout);
    
    return 0;
}

Что делает LinearFit()?
Эта функция принимает 2 аргумента:

1. Матрицу с данными, которые надо аппроксимировать.
   
2. Модельную функцию (точнее, функтор)

Матрица с данными имеет следующую структуру:

x[0,0],   x[0,1],   ..., x[0, M-1],   y[0]
x[1,0],   x[1,1],   ..., x[1, M-1],   y[1]
                    ...
x[N-1,0], x[N-1,1], ..., x[N-1, M-1], y[N-1]

где x — независимые переменные, а y — зависимая переменная.
Независимых переменных может быть несколько (M).

Модель представляет собой сумму линейно-независимых базисных функций с коэффициентами:

y = A1*f1(x1, x2, ..., xM) + A2*f2(x1, x2, ..., xM) + ... +  ADim*f2(x1, x2, ..., xM)

Функтор, используемый для задания модели должен иметь operator(), возвращающий значение некоторой базисной функции в некоторой точке. Он принимает 2 параметра:

1. номер базисной функции.
   
2. указатель на массив значений независимых переменных.

Кроме того, этот функтор должен иметь функцию Dim(), возвращающую размерность модели (кол-во базисных функций в ней).

Результатом применения МНК будут значения коэффициентов (A1, ... ADim).
Эти значения возвращаются из функции LinearFit() в виде матрицы [Dim * 1].