[iiv 17]

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

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

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

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

[Tanya 4]

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

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

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

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

[amaora 20]

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) Проблему большого набора данных этот метод не решает, все так же будет зависеть от реализации процедуры разложения.

 

 

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

[iiv 17]

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 хотя бы несколько раз обращаться, а не получать ее по строкам.

[amaora 20]

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

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

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

 

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

 

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

[iiv 17]

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 младшие собственные числа и не мудрствовать лукаво - я вам это много раз советовал, но вы, как я понимаю, так и не попробовали, а зря.

[amaora 20]

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]).

[iiv 17]

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

 

[amaora 20]

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

Так вот, а сейчас хочется решать системы высокой размерности для которых существует много решений не сильно отличающихся по минимизируемой квадратичной метрике. И пока не очень понимаю как например выбирать параметр в регуляризации Тихонова. А может быть надо минимизировать норму искомого неизвестного вектора b, есть такая функция в lsqminnorm в matlab, но мне нужна своя реализация. Здесь проблема, что непонятно какой критерий "правильного" решения выбрать. Хотелось бы, чтобы метод не был сильно сложнее обычного МНК через QR разложение.

Задачи теперь из области машинного обучения.

On 3/30/2021 at 11:10 PM, iiv said:

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

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

А для случая QR разложения как это будет выглядеть?
 

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

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

b = Rx \ S

трансформируется в,

[Q, R] = qr([[X Z]; eye(nx,nx)*la; zeros(nx,nz)])

// где, R = [Rx S; 0 Rz]
//      la = epsilon * max(abs(diag(Rx)))

b = Rx \ S

так?

[iiv 17]

4 hours ago, amaora said:

А для случая QR разложения как это будет выглядеть?

я боюсь, что до конца не понял нотацию матлаба, поэтому запишу обычными формулами: у вас есть без регуляризации \(b = R_x^{-1} s\), теперь это будет в виде \( b = (R_x^T R_x + \epsilon I)^{-1} R_x^T s \)

По поводу выбора регуляризационного параметра. Если вы взяли решение \( b_0 \) при нуле и при очень маленьком значении эпсилон этого параметра \( b_{\epsilon} \) и при удвоенном значении эпсилон этого параметра \( b_{2\epsilon} \) , и сравнили по норме и \( || b_{2\epsilon} - b_0 ||_2 > 2 || b_{\epsilon} - b_0 ||_2 \), тогда регуляризатор нужен, и обычно его оптимум достигается при таком эпсилон, когда это неравенство перестает быть верным, хотя и иногда немного раньше. Как угадать это эпсилон без перебора не дорого - готового решения нет. Если Вы можете для вашей \( R_x \) вычислить сингулярное разложение, то тогда такой поиск будет Вам стоить \( O(N^3 + kN)\), где N размерность \( R_x \), а k - двоичный логарифм машинной точности, если не можете, то \( O(kN^3)\).

[amaora 20]

Зачем делать квадрат исходной матрицы? Не знаю как здесь формулы добавлять, не буду пока отвечать.

Изменено 2 февраля, 2023 пользователем amaora

[iiv 17]

5 hours ago, amaora said:

Зачем делать квадрат исходной матрицы? Не знаю как здесь формулы добавлять, не буду пока отвечать.

Что вам запрещает сделать сингулярное вашей Rx, тогда в моих формулах вместо Rx будет стоять диагональная матрица D, и у вас будет полный контроль точности при вычислении этих формул: \( b = V (D + \epsilon D^{-1})^{-1} U^T s \) и никаких ошибок у вас не набежит. Без сингулярного только квадрат как я написал, или ранк ревеалинг с отбросами векторов на основе критерия, схожего с регуляризацией по Тихонову, но, боюсь, что с последним будет очень много танцев с бубнами.

[amaora 20]

21 hours ago, iiv said:

Что вам запрещает сделать сингулярное вашей Rx

Рассматриваю все варианты, вычислительная эффективность имеет значение. Вот так я добавляю новые уравнения \( x b = z \) в систему,

\( R' = qr\left ( \begin{bmatrix}
R_x && S_{xz} \\
0 && R_z \\ 
x && z
\end{bmatrix} \right )
,
R = \begin{bmatrix}
R_x && S_{xz} \\
0 && R_z
\end{bmatrix} \)

после накопления нахожу решение,

\( b = R_x^{-1} S_{xz} \) .

А теперь хотелось добавить в систему уравнений регуляризующие строки,

\( R' = qr\left ( \begin{bmatrix}
R_x && S_{xz} \\
0 && R_z \\ 
\epsilon I_x && 0_z
\end{bmatrix} \right ) \) .

Но будет ли это эквивалентно регуляризации Тихонова.

 

 

[amaora 20]

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

\(R(k+1) = qr(R(k)^T,0)\).

Повторяем транспонирование и QR факторизацию несколько раз, смотрим значения R(1,1) и R(n,n), там лежат оценки максимального и минимального сингулярных чисел. Или можно брать l2 норму от первой строки и последнего столбца, а сходимость отслеживать по тому насколько обнулялись недиагональные элементы. По тестам со случайными матрицами получается 4-5 итераций мне достаточно.

Остался только вопрос, действительно ли такой процесс будет всегда сходиться? Не находил нигде описания именно в таком виде, у меня упрощённый процесс без Q и без сложных поворотов с двух сторон. А в этих методах SVD часто бывает плохая сходимость в некоторых случаях. Но обычно все пытаются найти полный набор сингулярных чисел, мне же достаточно экстремальных.

Изменено 18 февраля, 2023 пользователем amaora

[amaora 20]

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

1. Находит оптимальное решение по критерию минимизации L2-нормы невязки (наименьшие квадраты);

2. Использует QR разложение реализованное с помощью быстрых вращений Гивенса (без вычисления квадратного корня);

3. Использует компенсацию потери точности из-за округления с помощью каскадного обновления матрицы R (увеличивает требуемое количество памяти, почти ничего не стоит по производительности);

4. Не требует хранения всех входных данных, можно отдавать по одной строке;

5. Есть оценка СКО невязки и min/max сингулярных чисел;

6. Есть L2-регуляризация для систем с недостатком ранга (но это можно делать и простым добавлением специально выбранных уравнений);

7. Простая реализация на C, из внешних зависимостей функции fabsf и sqrtf;

 

Проблема из п.3 проявляется на больших наборах данных, для single точности можно попробовать >1M уравнений, например в octave:

octave:17> X = single(randn(5000000,5));
octave:18> Y = single(X * [1; 1; 1; 1; 1]);
octave:19> ols(Y,X)
ans =

   1.0057
   1.0057
   1.0058
   1.0058
   1.0057

 

Ссылка на код - https://github.com/rombrew/lse