Имеется матрица комплексных значений размером 4 на 4.
Необходим код на Си решающий задачу определения собственных векторов указанной матрицы.
Дополнительные сведения:
1. Собственные числа матрицы уже определены с помощью алгоритма QR.
2. Система построена на сигнальном процессоре Аналог девайс DSP TigerSHARC.

Помогите ссылкой на удобный в реализации на процессоре алгоритм или готовый код программы.

Задача получения собственных векторов к уже имеющимся собственным значениям возникает довольно часто. Особенно, когда нужны не все собвектора, а только для малого числа выбранных собзначений (обычно либо самых больших, либо самых маленьких). Опять же хочется сэкомить вычисления, рассчитывая на то, что определить избранные собвектора будет быстрее, чем их полный набор.

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

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

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

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

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

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

Кстати, так как не сказано, что матрица - эрмитова, у Вас есть гарантии, что там нет жордановых форм, или их Вам тоже надобно найти?

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

Это корреляционная матрица, симметричная, обращаемая, с комплексными значениями. Математическое моделирование показало, что QR алгоритм работает, правда сходится медленно. На Си написать можно, но если уже имеется готовое решение на Си то лучше им и воспользоваться, а не изобретать велосипед. Если вы в теме то не могли бы дать ссылки на Lapack в сорсах, f2c. метод Арнольди и пояснить как это скачать.

да, тогда она у Вас не симметричная, а эрмитова, жордановых форм у Вас не будет, все собственные значения должны быть не отрицательны, и Арнольди Вам попросту будет не нужен.

Самое простое, скачать, скомпилить и радоваться http://www.netlib.org/lapack может Вам повезет и Вы найдете уже готовую сборку лапака для TigerSHARC, на раз гуглом я не нашел, но об этом много пишут, то есть должно быть.

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

При Вашей маленькой размерности может подойти и быть почти оптимальным по скорости еще один вариант:

Пусть Ваша матрица A=A_1,
выполним для нее последовательно следующие операции
A_{i+1} = A_i * A_i / ||A||_2^2
||A_i||_2^2 - сумма квадратов всех матричных элементов.

За число итераций, не более машинной точности (53 для двойной точности) Вы сойдетесь к одноранговой матрице, любая строка которой после нормировки равна собственному вектору, соответствующему максимальному собственному значению.

Отнормировав собственный вектор v_1 и посчитав собственное его значение по формуле l_1=v_1^* A v_1 Вы можете вычесть из исходной матрицы l_1 v_1 v_1^* и найти следующий, второй собственный вектор. Также дальше можно поступить и с третьим вектором, а для четвертого Вам и в квадрат-то возводить не надо будет.

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

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

Во-первых, что такое у вас ||A||_2? Я так понимаю, вы хотели отнормировать матрицу, но тогда, вообще-то, надо указывать какая норма матрицы считается.

Во-вторых, ваша итерационная процедура, если я её правильно понял, сойдётся в конце концов к диагональной матрице с собственными числами на диагонали. А причём тут собственные вектора?

В-третьих, какой lapack, какой blas?! Вы вообще когда-нибудь програмили ЦПОСы для real-tim'а? Открою вам небольшой печальный секрет: на ЦПОС можно запрограмить "примерно в 100 строк за вечер" только то, что использует прилагаемую разработчиком математическую (или сигнальную) библиотеку. Ибо, всё иное либо: а) приведёт к совершенно не realному времени (а нахрена тогда нужен DSP'шник - считайте на калькуляторе), б) наём дорогостоящего программера, который вашу функцию будет ногами и руками упихивать в прокрустово ложе конкретной архитектуры на асме.

2 Kluwert, не буду с Вами дискутировать, видя у Вас полное отсутствие знаний по вычислительной математике, тому пример, Ваше Высказывание:


Без уважения к Kluwert и с уважением к остальным

ИИВ

Спасибо за советы. Я так понял что процесс A_{i+1} = A_i * A_i / ||A||_2^2, где ||A_i||_2^2 - сумма квадратов всех матричных элементов должен быть многократно повторен чем больше повторений тем точнее результат как и при QR алгоритме. Здесь учитывается что матрица комплексная? Потому что алгоритм Грама-Шмидта в QR алгоритме отказался работать с комплексными числами. Здесь пришлось комплексную матрицу преобразовать в матрицу действительных чисел но уже в два раза большего размера и по завершении QR алгоритма результат преобразовать в комплексную матрицу. Насчет сравнительной эффективности этих методов вот в чем вопрос.



Спасибо за советы. Я так понял что процесс A_{i+1} = A_i * A_i / ||A||_2^2, где ||A_i||_2^2 - сумма квадратов всех матричных элементов должен быть многократно повторен чем больше повторений тем точнее результат как и при QR алгоритме. Здесь учитывается что матрица комплексная? Потому что алгоритм Грама-Шмидта в QR алгоритме отказался работать с комплексными числами. Здесь пришлось комплексную матрицу преобразовать в матрицу действительных чисел но уже в два раза большего размера и по завершении QR алгоритма результат преобразовать в комплексную матрицу. Насчет сравнительной эффективности этих методов вот в чем вопрос.

не совсем так, QR не обязан сойтись за фиксированное число итераций, а тот алгоритм, что я предложил сойдется (и точнее Вы все равно ничего не получите) не более, чем за 53 итерации для двойной точности, и не более 23 итераций для одинарной точности.



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


первый раз такое поведение Грамма-Шмидта вижу, программировал его с и без реортогоналицации сотни раз для всевозможных матриц, в основном для комплексных - все было хорошо. А Вы случаем, комплексное сопряжение не забывали где надо ставить? А, если честно, городить огород с 4х4 матрицей через Грамма-Шмидта - это перебор - Вам достаточно запрограммировать Гивенса, трижды его применить - получить трехдиагональную матрицу, а потом исчерпывать второе поддиагональное значение тем же Гивенсом. Если оное очень велико, то вначале матрицу переставить, снова применить три раза Гивенса, и гарантированно исчерпать полученное второе поддиагональное. Полученные две подматрицы размера 2х2 решить в лоб через квадратное уравнение, восстановить все необходимые вращения назад и радоваться длинному, но быстрому алгоритму. Или все-таки не полениться и скомпилить лапак и не ловить месяцами свои баги и мучиться с не понятно почему плохой сходимостью. Люди на STM32F103 GSL с лапаком компилили, а Вы боитесь.

2 Xenia

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

Удалось?
Ай да Бересфорд... А за пол-итерации можно?

Все удалось но только ж...пу отсидел. И еще выполнение операций в фомате float не позволяло пеленговать сигналы приходящие с малым угловым разрешением. И только переход на long double дал результат.

Можно было попробовать в двойной, но с использованием точных методов обращения матриц, или методами удвоения точности при суммировании и умножении (Кнут, Кахан, Бабушка и др.).
По теме, есть книга: Бибердорф Э.А., Попова Н.И. Гарантированная точность современных алгоритмов линейной алгебры.
Книга по работам авторов и материалам Годунова. Описаны алгоритмы вычисления векторов с гарантированной точностью. В том смысле, что кроме результата есть оценка точности вычисления этого результата.

Оно и без программирования было ясно, что на float'е в матричной алгебре далеко не уедешь. Но процессор у вас прежний остался, TigerSHARC? Какой разрядности у него doublе и long double? Аппаратная арифметика или эмуляция?

P.S. Давно (в 2008 году) у нас на форуме тема была, где Sharc-процессоры сильно ругали по поводу плавучки:
http://electronix.ru/forum/index.php?showtopic=55332
Насколько это справедливо к вашему TigerSHAC? Что-то сдинулсь с тех времен?
Может быть вам на другой процессор перейти, где double64 аппаратно поддерживается?

Это DSP и арифметика у него аппаратная. Платформу TS201 выбирали другие специалисты и кстати не прогадали если
учесть что в следующем году Миландер по программе импортозамещения выпускает процессор 1967ВЦ2Т и обещает
обратную совместимость с ADSP-TS201 как по исполняемому кода, так и по работе со средой разработки VisualDSP,
Минимальное изменение системной платы при миграции с TS201 на 1967ВЦ2Т.
Использовался Си компилятор поставляемый со средой VisualDSP и он соответсвует стандартну. Со средой разработки
поставляется библиотека Синапс для работы с комплексными числами complex_long_double и complex_float, матрицами, БПФ,
цифровые фильтры и множество разных функций. Но нужно вниметельно смотреть на версию среды разработки так как
несмотря на то что названия функций в зависимости от арифметики оканчиваются либо на f или d в старых версиях это
всегда f.

Ничего, новая отрастёт... А вообще, это прогресс.

У вас эрмит 4х4, можно вас попросить для примера привести по паре матриц для больших и малых углов.

А разве в TS201 long double отличается от double?