Вычисление собственных векторов Опции

[mihalevski]

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

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

[iiv]

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

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

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

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

[mihalevski]

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

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

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

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

[iiv]

Это корреляционная матрица, симметричная, обращаемая, с комплексными значениями.

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

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

[mihalevski]

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

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

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

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

[iiv]

должен быть многократно повторен чем больше повторений тем точнее результат как и при QR алгоритме.

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

Здесь учитывается что матрица комплексная?

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

Потому что алгоритм Грама-Шмидта в QR алгоритме отказался работать с комплексными числами.

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

2 Xenia

Об этом еще Парлетт писал в своей книжке "Симметричная проблема собственных значений". И, кажется, именно он это свойство впервые обнаружил.

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