Jump to content

    
Sign in to follow this  
amaora

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

Recommended Posts

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

похоже ни кого таких задач нет.

второй семестр нумерики любого ВУЗа по специальности "вычислительная математика" причём пройти на следующий уровень нельзя, если не умеешь такие задачи решать.  Поэтому я реально не понимаю ТС, что он всячески уповает только на точность вычисления матриц, забывая что он возводит в квадрат обусловленность, пренебрегает оценкой разброса сингулярных чисел и согласованности правой части с левыми сингулярными векторами.

В теме много ключевых слов вам насоветовали, но вы к ним не прислушиваетесь, а жаль, и, именно из-за этого вам и кажется, что вы "переоткрываете Америку".

Share this post


Link to post
Share on other sites

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

Вот можно под таким углом посмотреть... 

Будем её решать генетическим алгоритмом... Параметры от времени как будут меняться? 

Предлагаю познакомиться с такими словами - бифуркациями, переход к хаосу, аттрактор, биллиард Синая.

Share this post


Link to post
Share on other sites
On 3/16/2021 at 4:04 PM, iiv said:

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

Видимо о чем то разном говорим. Исходная постановка задачи, есть уравнения:

Quote

x(i) * b = z(i)

количество измерений i много превышает размерности векторов в уравнении. Переопределенная система, которую надо решить по критерию минимизации суммы квадратов невязок. Пусть это будет называться частным случаем регуляризации Тихонова, если так хочется.

Quote

sum_i || x(i) * b - z(i) ||^2 -> min

Известно "простое" решение (1), где X и Z матрицы составленные из исходных векторов-строк x(i) и z(i).

Quote

b = inv(X' * X) * X' * Z

Известно решение (2) с помощью QR разложения, блочный метод (наиболее удобная формулировка как мне кажется):

Quote

 

[Q, R] = qr([X Z])

// где R = [Rx S; 0 Rz]

b = S \ Rx

 

В смысле потери точности при плохо обусловленной матрице X, очевидно, что выигрывает метод (2) ортогональных преобразований. Но это не главная проблема из-за которой я начал тему. Оба метода деградируют при большой размерности набора данных при любой обусловленности матрицы X. Многое зависит от реализации процедуры QR разложения в методе (2) или метода накопления матрицы X' * X в методе (1).

* дополнительное условие задачи: нет доступа ко всем строкам X и Z одновременно, данные поступают последовательно, хранить их все негде.

 

Анализ через SVD меня не заинтересовал, по следующим причинам:

1) Редко встречаются плохо обусловленных задачи которые нужно решать в том виде как есть, и нет возможности набрать более информативные данные;

2) В связи с дополнительным условием мне нужна будет процедура обновления разложения при поступлении нового измерения, как я понял это не слишком просто;

3) Проблему большого набора данных этот метод не решает, все так же будет зависеть от реализации процедуры разложения.

 

 

Edited by amaora

Share this post


Link to post
Share on other sites
On 3/17/2021 at 6:10 PM, amaora said:

Видимо о чем то разном говорим. Исходная постановка задачи, есть уравнения:

количество измерений i много превышает размерности векторов в уравнении. Переопределенная система, которую надо решить по критерию минимизации суммы квадратов невязок.

Ага, верно. Вроде о том же. Я про сингулярные числа говорю только из-за того, чтобы понять какие они у вас. Я не советовал их вычислять.

 

 

Тут конкретика двояка (даже трояка), но, для начала условимся:

sum_i || x(i) * b - z(i) ||^2 -> min

и рассмотрим A = [X Z]

 

 

Если у этой матрицы разброс сингулярных значений меньше обратного из корня из машинной точности, то надо решать b = inv(X' * X) * X' * Z  и не морочить себе голову.

 

Если у этой матрицы разброс сингулярных значений больше обратного из машинной точности, то надо выбрасывать часть сингулярных чисел, иначе вы хоть даже в точной арифметике будете получать погоду на марсе вместо решения. Для этого - применимы методы регуляризации по Тихонову, и классикой является, например, такое решение b = inv(X' * X + a I) * X' * Z, где I - единичная диагональная, а скаляр a - величина близкая к отбрасываемым сингулярным значениям. Я не призываю сингулярные числа искать, надо иметь их оценку, чтобы регуляризовать.

 

Если же есть желание решать в промежуточном диапазоне, то тут только метод построения X' * X в удвоенной точности, и потом обращение с учетом того, что числа даны в удвоенной точности. Идеальным вариантом тут подходит

1. вначале построение в удвоенной точности A  = X' * X, и B = X' * Z,

2. вычисление Холецкого с выбором диагонального элемента или блока A = L L'  (тоже в удвоенной точности),

3. вычисление сингулярного разложения L = UDV' уже в обычной точности,

4. перезапись A = U D^2 U' и использование его далее, то есть умножение в удвоенной точности C = U' B, вычисление F = inv(D^2) C, и окончательное умножение в обычной точности b = U F.

 

Все остальное может быть применимо только если вы можете к самой матрице X хотя бы несколько раз обращаться, а не получать ее по строкам.

Share this post


Link to post
Share on other sites
On 3/29/2021 at 12:32 AM, iiv said:

Все остальное может быть применимо только если вы можете к самой матрице X хотя бы несколько раз обращаться, а не получать ее по строкам.

Сделать какие-то оценки чтобы понять, что получился мусор а не решение было бы полезно, но пока обойдусь. Априорно предполагаю хорошую обусловленность системы, это обеспечивается в большинстве случаев, т.к. я сам задаю "зондирующий сигнал X". Ошибка копится при агрегировании большого количества данных в одну матрицу, а делать это в двойной точности слишком дорого.

 

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

 

https://sourceforge.net/p/liblse/code/ci/default/tree/README.md

Share this post


Link to post
Share on other sites
34 minutes ago, amaora said:

Ошибка копится при агрегировании большого количества данных в одну матрицу, а делать это в двойной точности слишком дорого.

Вы когда умножаете, не держите двойную точность, а используйте ее только при сложении. У умножения ошибка - последний бит, а при сложении всяко может быть.

for(i=0; i<N; i++)

{ T = x x^T;  //  x и T - в одинарной точности

   A +=T; // A - в двойной точности

}

Если двойной точности нет, то пользуйте сумму двух float, сдвинутых друг от друга на машинную точность.

float s1, s2, a; // s1 больше s2 на обратную величину мантиссы. К этой паре (s1,s2) мы хотим добавить "a", для этого выполнив:

float t1=a1+a; // здесь t1,t2,t3 - промежуточные временные переменные

float t2=a-(t1-a);

t2+=s2;

s1=t1+t2;

float t3=s1-t1;

s2=(t1-(s1-t3))+(t2-t3); // по окончании s1,s2 содержит сумму a и того, что было до этого в s1,s2

 

А лучше просто выбрасывать у X' X младшие собственные числа и не мудрствовать лукаво - я вам это много раз советовал, но вы, как я понимаю, так и не попробовали, а зря.

Share this post


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

Если двойной точности нет, то пользуйте сумму двух float, сдвинутых друг от друга на машинную точность.

Раньше так и делал, выше приводил тест такого метода на octave. Ошибка агрегации хорошо компенсируется. Но это все же квадраты от исходных данных, и деградация точности тоже по квадрату от числа обусловленности. Смысл в этом методе только как в наиболее вычислительно простом, если делать LDL разложение, то можно и без квадратных корней обойтись.

 

1 hour ago, iiv said:

А лучше просто выбрасывать у X' X младшие собственные числа и не мудрствовать лукаво - я вам это много раз советовал, но вы, как я понимаю, так и не попробовали, а зря.

Про регуляризацию с параметром не понял как это надо использовать. Предлагаете под каждую задачу анализ проводить? Так и приедем к тому чтобы SVD делать.

 

On 3/29/2021 at 12:32 AM, iiv said:

Если у этой матрицы разброс сингулярных значений меньше обратного из корня из машинной точности, то надо решать b = inv(X' * X) * X' * Z  и не морочить себе голову.

А я анализирую число обусловленности только X матрицы а не [X Z], так даже непонятен смысл, у меня всегда огромные числа cond([X Z]).

Share this post


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

Про регуляризацию с параметром не понял как это надо использовать. Предлагаете под каждую задачу анализ проводить? Так и приедем к тому чтобы SVD делать.

не, тут все просто.

Если вам заранее известно, что исходные данные имеют погрешность epsilon, то в большинстве случаев достаточно трансформировать b = inv(X' * X) * X' * Z в

b = inv(X' * X + epsilon^2 * max_L * I) * X' * Z, где max_L - максимальное собственное значение X' * X.

8 minutes ago, amaora said:

А я анализирую число обусловленности только X матрицы а не [X Z].

а вот это - зря. Во-первых, надо сравнивать обусловленность X матрицы с некоторой абстрактной обусловленностью [X Z] без самого маленького сингулярного числа.

По хорошему [X Z] - должна быть вырожденной, с единственным нулевым сингулярным числом. Вот есть его выбросить, и выбросить все сингулярные числа меньше погрешностей во входных значениях, у вас получится очень интересная матрица, она вроде и вырожденная, но у нее есть решение с минимальной нормой.

Чтобы не мудрствовать и не считать SVD (его всегда неудобно считать) как раз и играются со сдвигом inv(X' * X + epsilon^2 * max_L * I) * X' * Z - и, в большинстве случаев, имея даже половинную точность (как в алгоритмах распознавания и искусственного интеллекта) имеют приемлемое решение.

 

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