Вопрос: есть ли в GSL (ну, пусть даже не GSL, а BLAS или даже LAPACK) методики решения таких уравнений?

Ответ: понятия не имею но я бы посоветовал набросать метод в python затем проверить и потом заменить тип на комплексные числа после чего еще раз проверить все и не морочить себе голову.

Умник! Ты в курсе, чем система с комплексными числами отличается от системы с действительными?

эм... так если искомый вектор x - вектор действительных чисел, то система первращается в A'x = b', где A' и b' - матрица и вектор действительных чисел соответственно.

Умник! Ты в курсе, чем система с комплексными числами отличается от системы с действительными?

ВНЕЗАПНО там обычные числа а там комплексные или не?

Вообще-то это даже справедливо и для комплексных x.

Система будет { Re(A)*Re(x) - Im(A)*Im(x) = Re(b), Re(A)*Im(x) + Im(A)*Re(x) = Im(b) (поэлементное взятие)

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

A — это базис векторов, B — искомое векторное поле (вектора двумерные, поэтому я и решил их комплексными числами заменить).

Естественно, коэффициенты x должны быть действительными (т.к. по сути мы ищем оптимальные коэффициенты разложения поля B по базису A).

И никак не могу понять, как мне это считать!

Да, получаем две системы, но в них одни и те же коэффициенты X.

И не надо умиляться! Как мне эти X состыковать? Если я по-отдельности решу для действительной и мнимой части, получу два разных вектора, которые будут ни к селу, ни к городу!

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

Всё.

Ну так в этом случае у тебя получается система { Re(A)*x = Re(b), Im(A)*x = Im(b) Правда матрица уже не ортогональная. Не знаю, важно ли это.

Всё

А вот ничего подобного! В этом случае будет то же самое, что по-отдельности решать, т.к. длина ответа будет в 2 раза больше нужной.

В том-то и "прикол", что хрен его знает, как это решить.

Эдди, ты совсем плохой стал.

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

Возьми {{1+i, 2+i}, {1-i}, {2-i}} и честно распиши итоговую 4x4 матрицу. Она _не_ будет блочной, её надо будет приводить к такому виду.

А про длину ответа - вообще детский сад. Про инвариантность линейных пространств забыл?

И если вообще матрицы прямоугольные, о каком «решении» ты вообще говоришь? Если у неё строк меньше чем столбцов - то у тебя не базис, а говно, если же наоборот, то базис какого-то подпространства, и для b вне его решений у твоей задачи не будет вовсе. Можно искать решение для проекции b на него, задача будет другая, но решаться будет ровно также.

мало данных (ошибка зверская вылезает)

Что это значит? Размерности базиса не хватает?

Ты не можешь решать по-отдельности

Я тебе о чем толкую? Ясен пень, по-отдельности нельзя. А ты именно по-отдельности и предлагаешь!

о каком «решении» ты вообще говоришь?

Ты про наименьшие квадраты слышал? Которые невязку к минимуму устремляют.

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

При равномерном распределении и количестве точек где-то эдак хотя бы около 10000 (100х100) получается очень даже прелестно обычным разложением по базису восстановить (точность восстановления коэффициентов получается около 0.1%, точность восстановления данных — тысячные доли процента), в этом случае МНК дает намного худшие результаты.

А вот когда данные распределены так, как у меня, "нормальное" разложение дает сильные ошибки (из-за больших дырок), да и понятно, что всего-то 256 точек физически не позволят получить коэффициенты разложения с точностью выше нескольких процентов! Вот и хочу МНК.

Это не решение системы линейных уравнений. Это решение ровно того, о чём я написал. GMRES тут не при чём.

Блин, ну ты что, никогда не решал системы методом наименьших квадратов?

Скажем, есть у тебя тысяча измерений какой-то величины, ты знаешь примерную зависимость и при помощи МНК (в октаве это оператор \) считаешь.

Для разреженных матриц GMRES - классика, хотя сейчас BCG с толковым предобуславливателем лучше.

А я привык работать с формальными определениями задачи и методов её решения, а не твоё карго знание операторов Octave.

Ну смотри: разложить, как ты предлагаешь, нельзя — вместо ответа получим говно. Всякие классические решаторы (QR, LU, хаусхолдер и т.п.) либо не работают в GSL с комплексными числами, либо требуют квадратную невырожденную матрицу. Я маленько нагуглил, что в LAPACK что-то подобное должно быть, но глянул синтаксис этой библиотеки и быстренько жамкнул ctrl+w: это вообще жесть жуткая, фортран какой-то поверх сей!

Для начала у нас вместо задачи говно. Если рассмотреть чистую линейную алгебру - у тебя Ax=b (где всё уже действительное), с параметризацией: x_{N+i} = 0, так?

фортран какой-то поверх сей!

А ты быстрый.

http://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)#The_general_p...

Последняя формула, не?

Если рассмотреть чистую линейную алгебру, то у меня прямоугольная матрица A высотой Sz (количество точек) и шириной pmax (количество аппроксимирующих полиномов, составляющих ортонормированный базис), которая умножается на вектор X (коэффициенты перед полиномами), давая мне векторное поле Y: AX=Y. Каждый столбец матрицы A — это значение i-й базисной функции в точках, где определено мое векторное поле Y. Понятно, что столбец тоже является векторным полем, т.е. комплексный. При этом коэффициенты X — действительные.

Методом наименьших квадратов мне нужно найти X, которое дает AX=Y+d, где d — комплексная невязка, а X такое, что сумма квадратов модулей d минимальна!

Что уж тут непонятного?

У тебя же эти «точки» — координаты исходного вектора b в некотором пространстве (в котором также определён базис A). Как ты берёшь их произвольное число? У твоего пространства размерность какая?

Если ты рассматриваешь исходный вектор b как проекцию «истинного вектора» на область пространства, определяемую известными координатами (т.е. некоторые координаты известны, остальные берём равными нулю) то как от него можно перейти обратно к «истинному» вектору? Если координаты независимы, то у «истинного» вектора неизвестные координаты могут быть какими угодно. Они же по отношению к известным координатам ортогональны в самом что ни на есть прямом смысле этого слова.

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

тогда в чем проблема?

A^T A x = A^T b => x = A^T b.

То, что ты привел — обычное разложение функции по базису. А у меня точек слишком мало и расположены они неравномерно. Я же уже говорил.

x = A^T b

Почти, x = (A^T A)^-1 b

Ты получил нелинейную задачу по минимизации квадратов сумм вида |\sum a_{ij} x_i - c_i| + |\sum b_{ij} x_i - d_i|, где A = Re(A), B = Im(A), c = Re(b), d = Im(b).

Так что хрен тебе, а не linear least squares, решать будешь общую задачу.

x = A^T b

Почти, x = (A^T A)^-1 b

перечитай еще раз. A - ортогональна.

Методом наименьших квадратов мне нужно найти X, которое дает AX=Y+d, где d — комплексная невязка, а X такое, что сумма квадратов модулей d минимальна!

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

Тогда проблема не в том, что X плохо вяжется с Y, а в том, что ты ищешь решение для X и Y', а Y' плохо вяжется с Y.

Как ты берёшь их произвольное число?

Поясню: это я решил сегодня добить-таки своих несчастных Церников. Точки — координаты пятен на изображении, вектора — направление лучей света, образовавших эти пятна (т.е. нормали к волновому фронту). Мне нужно восстановить волновой фронт. Матрица A — полиномы Жао в координатах пятен.

Раскладывая векторное поле (нормали к волновому фронту) по векторным полиномам Жао, я получаю путем несложных преобразований коэффициенты Цернике, которые позволяют мне восстановить волновой фронт. Вот и все. Элементарщина, но хрен посчитаешь.

То, что ты привел — обычное разложение функции по базису. А у меня точек слишком мало и расположены они неравномерно. Я же уже говорил.

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

А пока мне неясно где проблема. Ведь транспонирование и умножение матрицы на вектор не дают дополнительных погрешностей.

☹

Ну неужто низзя?

Получается, что базис у тебя таки неполный

Естественно: я же ограничиваюсь первыми pmax полиномами! Не буду же я бесконечное количество полиномов фигачить!

неважно полный, или неполный. Проекцию на это пространство получаешь обычным умножением транспонированной матрицы.

Базис неполный (полный — бесконечный).

Представь, что мы берем функции sin(x), sin(2x) и sin(3x) и пытаемся ими аппроксимировать функцию, заданную точками с x=0, x=0.1, x=pi/2-0.1 и x=pi/2+0.1.

Дык, я уже со скалярными функциями проверял: МНК дал мне невязку в 10^{-23} на подопытном базисе, а непосредственное разложение — 10^{-3}. Как тебе разница?

где погрешность? В векторе y(x) - аппроксимационных точках?

Дык, я уже со скалярными функциями проверял: МНК дал мне невязку в 10^{-23} на подопытном базисе, а непосредственное разложение — 10^{-3}. Как тебе разница?

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

А, забыл добавить: в точках же не точное значение векторного поля, а с некоторой погрешностью.

Думаю, что в таком случае полигоны Жао твои — полное говно. Точно помню, что отбрасывание верхушки ряда в нормальных базисах даёт проекцию, которая наиболее близка к исходной сумме ряда (раскладываемой функции) как раз в среднеквадратическом смысле, будь то разложение в линейном пространстве, будь то разложение по базису Фурье.

A - ортогональна

Чо?

Решение наименьшими квадратами прямоугольных систем с комплексными числами (комментарий)

У него просто прямоугольная матрица.

Ты понимаешь, что на неравномерном наборе точек по сути получится взвешенное разложение, где вес будет тем больше, чем выше плотность точек?

Если рассмотреть чистую линейную алгебру, то у меня прямоугольная матрица A высотой Sz (количество точек) и шириной pmax (количество аппроксимирующих полиномов, составляющих ортонормированный базис),

ортонормированный базис

У тебя полиномы нормальны или столбцы матрицы?

У него просто прямоугольная матрица.

из столбцов ортонональной матрциы A => A^T A = I.

Я выше говорил: каждый столбец матрицы A — это значения очередного полинома из базиса. Базис ортонормированный. Ортогональность доказана вышеупомянутым Жао, а нормирование выполняется на стадии вычисления значений этих полиномов в заданных точках.

Котелок уже закипел

В общем, я так понял, что никто не может предложить более-менее простой способ посчитать это дело при помощи МНК. Ладно, буду довольствоваться более хреновыми результатами.

Просто мне хотелось выпендриться и сделать точность на пару-тройку порядков выше, чем в коммерческой софтине (моя — под GPLv3).