[amaora 20]

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

 

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

 

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

Изменено 30 января, 2021 пользователем amaora

[amaora 20]

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

 

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

Изменено 6 марта, 2021 пользователем amaora

[Tanya 4]

У Вас линейная зависимость функции от параметров подгонки?

[amaora 20]

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

[Tanya 4]

15 часов назад, amaora сказал:

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

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

[amaora 20]

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

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

 

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

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

 

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

[iiv 17]

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.

[amaora 20]

У меня нет возможности хранить всю матрицу, надо аккумулировать измерения. Сделал в 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 не впечатляет.

Изменено 9 марта, 2021 пользователем amaora

[iiv 17]

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

2 hours ago, amaora said:

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

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

 

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

 

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

[amaora 20]

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

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

[iiv 17]

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 вам, к сожалению, не помогут.

[amaora 20]

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

 

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

 

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

 

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

 

Изменено 10 марта, 2021 пользователем amaora

[looser 8]

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

[amaora 20]

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

[iiv 17]

On 3/15/2021 at 1:43 PM, amaora said:

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

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