Почему Mldivide Matlab намного лучше, чем Dgels?

Question

Решить Ax = b. Настоящий двойник. A переопределено Mx2 с M >> 2. b равно Mx1. Я запустил кучу данных против mldivide, и результаты отличные. Я написал мексиканскую рутину с MKL LAPACKE_dgels, и это далеко не так хорошо. Результаты имеют массу шума, и основной сигнал едва там. Сначала я проверил процедуру по результатам примера MKL. Я искал документ mldivide (блок-схему) и вопросы SO. Все, что я нашел, это то, что Matlab использует QR-факторизацию для переопределенного прямоугольника.

Что мне следует попробовать дальше? Я использую неправильную процедуру LAPACK? Пожалуйста, помогите направить меня в правильном направлении.

Обновление: В пределах разницы с плавающей запятой E-15 по вектору решения Intel MKL LAPACKE_dgels имеет тот же результат, что и Matlab mldivide для реального двойного переопределения(прямоугольные) проблемы. Насколько я могу судить, это используемый метод QR.

Остерегайтесь остатков, возвращенных из этих dgels. Они не равняются b - топору. Многие из них близки к этому значению, а некоторые далеки от него.

Answer 1

Проблема была не в решении x, а в возвращенных остатках от DGELS. Выходные данные этой подпрограммы модифицируются на месте указателей входного массива. В MKL doc указано, что входной массив b перезаписывается выходным вектором x для первых N строк, а затем с остатками от N+1 до M. Я подтвердил это своим кодом.

Ошибка заключалась в выравнивании остатков b[N+1] с исходными входными данными b[1] и принятии дальнейших алгоритмических решений по этому вопросу. Правильное выравнивание остатка к исходному вводу составляет b[1] до b[1]. Первые N остатки недоступны;Вы должны вычислить их потом.

В документе не говорится, что они являются невязками как таковыми, а точнее:

остаточная сумма квадратов для решения в каждом столбце определяется каксумма квадратов модулей элементов от n+1 до m в этом столбце.