На самом деле, правильнее говорить не о методе, а о задаче наименьших квадратов, так как существуют несколько методов её решения. Сформулируем эту задачу. Пусть заданы матрица A размера n * m ( n >= m ) и вектор b размерности n. Требуется определить такой вектор x, чтобы выполнялось условие

    | Ax - b |2 = min    (1)

Обычно решение этой задачи сводят к решению системы, так называемых, нормальных уравнений:

    ATAx = ATb     (2)

Её можно получить взяв производную функции (1) по x, а затем приравнять её нулю. В этом случае мы получим квадратную (m*m) симметричную положительно определённую матрицу, которую можно решить, например, методом Гаусса или специальными методами для симметричных матриц. К сожалению, эта матрица часто имеет очень большое число обусловленности, что сказывается на точности результата. Поэтому лучше применять методы описанные ниже.

Эти методы ищут решение задачи (1) не приводя её к форме (2). Первый из них, основанный на ортогонализации Грама-Шмидта, представлен здесь в виде функции:

bool ortholin ( Matrix & a, CArrRef<double> b, ArrRef<double> x );

При этом матрица а изменяется ( в результате все столбцы матрицы будут ортогональными ).
Второй метод, с применением преобразований Хаусхолдера, может быть вызван при помощи функции:

bool lss_h ( Matrix & a, CArrRef<double> b, ArrRef<double> x );

Здесь также матрица а изменяется.
Вариант для комплекных чисел находится здесь.
В случае, когда надо решить несколько задач с одной матрицей, но разными векторами b, можно воспользоваться следующим классом:

class LSS_H
{
    Matrix & a;
    DynArray<nat> pivot;
    mutable DynArray<double> alpha;
// Запрет конструктора копии и оператора присваивания:
    LSS_H ( LSS_H & );
    void operator = ( LSS_H & );
public:
    LSS_H ( Matrix & a );
    bool solve ( CArrRef<double> b, ArrRef<double> x ) const; // b[a.nRow], x[a.nCol]
};

В серии моих испытаний второй алгоритм ( с применением преобразований Хаусхолдера ) гораздо чаще давал лучшую точность, чем первый ( ortholin ). Вероятно, это происходит потому, что реализация алгоритма с преобразованиями Хаусхолдера использует выбор ведущего столбца, а реализация алгоритма с ортогонализацией - нет. Впрочем, для большинства приложений оба варианта дают приемлимую точность. Эти же алгоритмы можно применять для случая, когда n = m, вместо метода Гаусса. При этом точность в сложных случаях будет лучше, правда, количество операций увеличится.

Эти алгоритмы были взяты мной из книги "Справочник алгоритмов на языке Алгол" ( Уилкинсон, Райнш ) и переписаны на С++ с небольшими изменениями. В частности я убрал процедуры уточнения результатов.

Описание класса Matrix находится здесь.
Описание шаблонов классов CArrRef, ArrRef и DynArray находится здесь.

Исходники находятся в файле mathem.cpp.