Jump to content

    
Sign in to follow this  
amaora

Численная реализация МНК

Recommended Posts

Часто возникает необходимость что-нибудь оценивать по измерениям. Для примера, расчёт импеданса по двум осям, расчёт момента инерции, нахождение калибровочных коэффициентов аналогового датчика угла. Все это нужно делать on-line на ... пусть cortex-m4f, с ограничением на время выполнения и занимаемую память. Сейчас использую примитивную реализацию с накоплением матриц-квадратов от входных сигналов $b = inv(X'*X) * X'*Y$. Вместо вычисления обратной использую LDLT разложение и прямую-обратную подстановку. Так же для получения достаточной точности начал при накоплении матриц $X'*X$ и $X'*Y$ использовать суммирование Кэхэна.

 

Сейчас подумал, может быть стоит перейти полностью на квадратно-корневую реализацию МНК. Хранить треугольную матрицу и делать операцию обновления (cholupdate) вместо суммирования. Не знаю нужно ли и можно ли в этом случае использовать суммирование Кэхэна. Хотелось бы получать методическую точность (для хорошо обусловленных задач) не зависящую от количества измерений и на уровне близком к машинному эпсилон. Размерности от 3 до 10-20.

 

Какие вы используете численные реализации МНК?

Edited by amaora

Share this post


Link to post
Share on other sites

Делал тесты в octave. Система из 3-5 неизвестных и ~15000 измерений. Квадратно-корневой метод через QR разложение быстрее теряет точность с ростом количества измерений, чем накопление с компенсацией Кэхэна и дальнейшее решение $b = inv(X'*X) * X'*Y$.

 

Вопрос в том, как применить какую-то компенсационную схему (подобно алгоритму Кэхэна) к QR разложению, для увеличения точности. Пока вижу возможность сделать аналог каскадного суммирования. Копим в промежуточной матрице порцию измерений, скажем 1000 измерений, затем этот кусок уже приведённый к треугольному виду "сбрасываем" в основную сумму. Так и сделаю наверно, размер порции только как-то адаптивно надо выбирать.

Edited by amaora

Share this post


Link to post
Share on other sites
15 часов назад, amaora сказал:

Да все задачи сводятся к виду $\b * x = y$, где \b искомый неизвестный вектор, а x и y известные сигналы.

Лучше бы тогда уж картинку.

Share this post


Link to post
Share on other sites

Вот здесь достаточно кратко про саму задачу, если я был непонятен.

https://ckrisgarrett.github.io/qr.html

 

Как я понял, то что мне нужно называется tall skinny QR, то есть QR декомпозиция "длинной матрицы" в которой гораздо больше строк чем столбцов. Но все, что мне удаётся найти по этому поводу направлено на параллельные вычисления, а не вопросы методической точности и компенсационных схем.

https://www.cs.purdue.edu/homes/dgleich/publications/Constantine 2011 - TSQR.pdf

 

Идея с разделением на квадратные блоки-матрицы и раздельной декомпозицией каждого куска понятна. Ожидаю, что точность будет деградировать аналогично каскадному суммированию. Но вот проблема, я использовал алгоритм Кэхэна из-за того ему не нужно держать в памяти весь массив данных а можно подавать измерения по одному, нужно лишь ещё немного памяти для остатков. Подобную схему для QR и хотелось бы найти.

Share this post


Link to post
Share on other sites
19 hours ago, amaora said:

Как я понял, то что мне нужно называется tall skinny QR, то есть QR декомпозиция "длинной матрицы" в которой гораздо больше строк чем столбцов. Но все, что мне удаётся найти по этому поводу направлено на параллельные вычисления, а не вопросы методической точности и компенсационных схем.

Чтобы немного больше понять про Вашу задачу, Вам надо вместо QR начать ее решать через SVD (вычислительная сложность больше только в несколько раз).

Пусть у Вас имеется $||Ax-b||_2^2$, где A содержит много меньше столбцов, чем строк. Сделайте для этой матрицы сингулярное разложение A=UDV^T, причем V - квадратное (в octave оно есть), и решение можно будет сформулировать как

x=VD^{-1}U^Tb,

Но тут важно следующее:

1. На диагонали D может быть что-то очень маленькое, но в соответствующей позиции вектора U^T b будет тоже что-то очень маленькое, в этом случае Вам повезло и этот "почти нуль" можно игнорировать.

2. На диагонали D может быть что-то очень маленькое, но в соответствующей позиции вектора U^T b будет что-то не маленькое, и из-за этого малые возмущения исходных данных будут приводить к огромным возмущениям результатов.

 

Если таки "почти нули" у Вас таки имеются (случай 2), то Вам надо добавить к D перед обращением единичную матрицу, умноженную на величину ошибки (eps) в исходных данных. В принципе это всегда правильнее делать.

 

Если писать в Ваших терминах через псевдообратную, (то есть Вы таки хотите использовать Холецкого или QR), то вместо x = inv(A^T A) A^T b Вам надобно считать  x = inv(A^T A + eps I) A^T b

 

PS: из-за того, что x = inv(A^T A) A^T b метод работает лучше, чем QR, у Вас, похоже, случай Nr 1.

Share this post


Link to post
Share on other sites

У меня нет возможности хранить всю матрицу, надо аккумулировать измерения. Сделал в octave тест точности разных методов:

 

1) Метод нормальных уравнений $b = inv(X'*X) * X'*Y$;

2) Метод нормальных уравнений но при накоплении матриц X и Y использовался алгоритм суммирования Кэхэна;

3) Обновление треугольной матрицы cholupdate (при этом содержимое матрицы аналогично QR разложению матрицы всех измерений [x y]);

4) Обновление треугольной матрицы cholupdate с использование каскадирования (три уровня).

5) QR разложение полной матрицы всех измерений;

6) SVD разложение полной матрицы всех измерений.

 

Spoiler

#!/usr/bin/octave -q

X = single(zeros(3,3));
Y = single(zeros(3,2));

Xk = single(zeros(3,3));
Yk = single(zeros(3,2));

X_rem = single(zeros(3,3));
Y_rem = single(zeros(3,2));

R = single(zeros(5,5));

Rk = single(zeros(5,5));
R1 = single(zeros(5,5));
R2 = single(zeros(5,5));
n1 = 0;
n2 = 0;

th0 = [1 2; 3 4; 5 6];

len = 300000;
Xfull = single(zeros(len,5));

for i=1:len

    % informativity signal
    xt = [1 1 1] + randn(1,3) * 0.0001;

    % measurement (with random noise)
    yerr = randn(1,2) * 0.0;
    yt = xt * th0 + yerr;

    x = single(xt);
    y = single(yt);

    Xfull(i,:) = [x y];

    % -- normal equations

    X = X + x' * x;
    Y = Y + x' * y;

    % -- normal equations (with Kahan summation)

    X_val = x' * x;
    X_fixed = X_val - X_rem;
    X_newsum = Xk + X_fixed;
    X_rem = (X_newsum - Xk) - X_fixed;
    Xk = X_newsum;

    Y_val = x' * y;
    Y_fixed = Y_val - Y_rem;
    Y_newsum = Yk + Y_fixed;
    Y_rem = (Y_newsum - Yk) - Y_fixed;
    Yk = Y_newsum;

    % -- QR

    R = cholupdate(R, [x y]');

    % -- QR (with cascaded accumulation)

    R1 = cholupdate(R1, [x y]');
    n1 = n1 + 1;

    if n1 >= 20
        for j=1:5
            R2 = cholupdate(R2, R1(j,:)');
        end
        R1 = single(zeros(5,5));
        n1 = 0;

        n2 = n2 + 1;

        if n2 >= 30
            for j=1:5
                Rk = cholupdate(Rk, R2(j,:)');
            end
            R2 = single(zeros(5,5));
            n2 = 0;
        end
    end
end

e_normal = norm(X \ Y - th0)
e_kahan = norm(Xk \ Yk - th0)
e_qr = norm(R(1:3,1:3) \ R(1:3,4:5) - th0)
e_qr_cascaded = norm(Rk(1:3,1:3) \ Rk(1:3,4:5) - th0)

[Q,R] = qr(Xfull,0);
e_qr2 = norm(R(1:3,1:3) \ R(1:3,4:5) - th0)

[U,S,V] = svd(Xfull(:,1:3),0);
e_svd = norm(V*inv(S)*U'*Xfull(:,4:5) - th0)
             

 

 

Результат для хорошо обусловленной системы.

Quote

e_normal =    7.3149e-04
e_kahan =    1.5232e-06
e_qr =  0.0032038
e_qr_cascaded =    6.5041e-06
e_qr2 =    5.3381e-04
e_svd =    6.2609e-04

И для плохо обусловленной.

Quote

warning: matrix singular to machine precision, rcond = 1.73611e-08
warning: called from
    lsq at line 86 column 10
e_normal =  4.0000
warning: matrix singular to machine precision
warning: called from
    lsq at line 87 column 9
e_kahan =  4.0000
e_qr =  0.97758
e_qr_cascaded =    1.9182e-04
e_qr2 =  19.624
e_svd =  26.617

 

Точность функций разложения в octave не впечатляет.

Edited by amaora

Share this post


Link to post
Share on other sites

Честно говоря, ничего не понял, что Вы написали :) но все равно попытаюсь прокомментировать.

2 hours ago, amaora said:

warning: matrix singular to machine precision, rcond = 1.73611e-08

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

 

Думаю, что мы тут все не экстрасенсы, чтобы угадывать, что конкретно означает ваши e_eps и остальные Ваши обозначения. Мне также не известно что конкретно Вы называете "алгоритмом суммирования Кэхэна". Напишите формулы, что Вы тут насчитали?

 

Вангую, что, скорей всего, ваши X по их физическому смыслу, линейно зависимы и, из-за этого, вы получаете такую чехорду с решениями и в каком-то из решений вы его случайно регуляризуете.

Share this post


Link to post
Share on other sites

Там полный код теста на octave раскрывается под надписью "Reveal hidden contents". Что не понятно?

Если не обнаружится ничего лучше, то наверно сделаю (или найду готовую) реализацию на C каскадного обновления QR, пригодную для использования на cortex-m4f.

Share this post


Link to post
Share on other sites
1 hour ago, amaora said:

Там полный код теста на octave раскрывается под надписью "Reveal hidden contents". Что не понятно?

только то, что конкретно вы там хотели посчитать и зачем.  Вы б математические формулы написали бы, что именно вы там считаете - всем стало бы понятно, не так ли?

 

В вашем случае поможет только то, чтобы считать inv(X^T X + alpha I) X^T y, и выбирать alpha в соответствии с тем, что я вам ранее посоветовал, не забывая, что у не отрицательно определенной матрицы X^T X собственные значения совпадают с квадратами сингулярных чисел X.

 

Ну и задуматься о том, что я тоже выше написал, по поводу того, что у вас в исходных X столбцы по физике похоже коллинеарны или хотя бы линейно зависимы со всеми вытекающими. То есть пока вы эту проблему не устраните, никакие QR да хоть с ранк-ревеалингом и в quad-double вам, к сожалению, не помогут.

Share this post


Link to post
Share on other sites

Есть исходные задачи, которые в физических терминах описаны. Там можно и нужно пытаться улучшить обусловленность при возможности, задавать наиболее информативный "зондирующий" сигнал.

 

А есть математический абстрактный решатель задачи наименьших квадратов. Здесь численный метод от которого требуется получать решение с наименьшей методической погрешностью. Я сравнил точности нескольких методов на числах одинарной точности, для уравнений с 6-ю неизвестными (матрица th0) по 300000 измерений. Два варианта входного "зондирующего" сигнала, информативный (столбцы в матрице A не коллинеарны) и слабо информативный (столбцы матрицы A отличаются в 5-м десятичном знаке). То, что метод нормальных уравнений развалился это ожидаемо, и я не планирую подставлять подпорки и пытаться это как-то исправить. Приемлемую точность выдал только метод каскадного обновления QR.

 

Сейчас у меня все задачи имеют хорошую обусловленность и реализован метод 2 (нормальные уравнения с компенсационным суммированием Кэхэна), точность получается приемлемая, это видно и по тесту. Но хотелось бы найти решатель на все случаи, в том числе и плохо обусловленные.

 

В случае плохой обусловленности будут сильно влиять погрешности измерений, и для нахождения точного решения надо делать много больше измерений. Здесь и возникает проблема деградации точности численного метода при большом количестве измерений. Найти по этой теме пока ничего не получилось.

 

Edited by amaora

Share this post


Link to post
Share on other sites

Просто мнк вам не подходит. Смотрите: в случае большого числа обусловленности при вычислении матрицы грама получаем квадрат числа обусловленности ну и последующие проблемы. В таких случаях используют приведение исходной матрицы к треугольной при помощи отражений (Хаусхолдер) или плоских вращений (Гивенс). Решение будет не в смысле мнк. Но вполне приемлемо.

Share this post


Link to post
Share on other sites

Это все МНК, при вычислениях неограниченной точности решение будет одно у всех методов. Уже набросал реализацию каскадного метода на C, для cholupdate использую вращения Гивенса. Всем спасибо, похоже ни кого таких задач нет.

Share this post


Link to post
Share on other sites
On 3/15/2021 at 1:43 PM, amaora said:

Это все МНК, при вычислениях неограниченной точности решение будет одно у всех методов

не верное утверждение, гуглите на методы регуляризации, в том числе Тихоновские, и читайте мои комментарии выше.

Share this post


Link to post
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

Sign in to follow this