Численный поиск обобщенных собственных векторов в Matlab

У меня есть матрица, такая как этот пример (мои фактические матрицы могут быть намного больше)

A = [-1 -2   -0.5;
      0  0.5  0;
      0  0   -1];

который имеет только два линейно независимых собственных значения (собственное значение -1 повторяется). Я хотел бы получить полную основу с обобщенными собственными векторами. Я знаю, как это сделать, используя функцию Matlab jordan в Symbolic Набор математических инструментов, но я бы предпочел что-то, предназначенное для числовых входных данных (действительно, с двумя выходными данными jordan не работает для больших матриц: «Ошибка в команде MuPAD: матрица сходства слишком велика».). Мне не нужна каноническая форма Жордана, которая, как известно, нестабильна в числовом контексте, мне нужна только матрица обобщенных собственных векторов. Существует ли функция или комбинация функций, которая автоматизирует это численно стабильным способом, или нужно использовать общий ручной метод (как стабильная такая процедура)?

ПРИМЕЧАНИЕ. Под «обобщенным собственным вектором» я подразумеваю ненулевой вектор, который можно использовать для дополнения неполного базиса так называемого дефектная матрица. Я не имею в виду собственные векторы, соответствующие собственным значениям, полученным в результате решения обобщенной проблемы собственных значений. используя eig или qz (хотя это последнее использование довольно распространено, я бы сказал, что его лучше избегать). Если кто-то не поправит меня, я не верю, что это одно и то же.


ОБНОВЛЕНИЕ 1 — пять месяцев спустя:

См. мой ответ здесь о том, как символически получить обобщенные собственные векторы для матриц размером более 82 на 82 (предел для моего теста матрица в этом вопросе).

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

ОБНОВЛЕНИЕ 2 — февраль 2015 г. Все вышеперечисленное остается верным, как было протестировано в R2014b.


person horchler    schedule 03.09.2013    source источник
comment
@natan: тогда я предлагаю вам пометить этот вопрос как дубликат.   -  person Jonas    schedule 04.09.2013
comment
@natan: Да, я видел этот вопрос/ответ. Это не дубликат. По сути, это ручная процедура, которая повторяет то, что показано во многих учебниках и в Википедии. Я сам знаю, как получить обобщенные собственные значения - я ищу что-то более автономное и использующее встроенные процедуры, чтобы быть осторожным с числовыми проблемами. Вы знаете что-нибудь или какие-нибудь алгоритмы? Кроме того, к вашему сведению, eigensys больше не существует в R2013a, необходимо использовать eig(sym(A)).   -  person horchler    schedule 04.09.2013
comment
@DavidHeffernan: ошибка возникает ровно 83 на 83 для конкретной тестовой матрицы.   -  person horchler    schedule 04.09.2013
comment
@DavidHeffernan: Что касается конечного размера, который меня интересует, эти матрицы представляют поведение узлов динамической нейронной сети, поэтому система может легко иметь более 100 измерений. Эти матрицы, якобианы, очень разрежены, поэтому было бы неплохо иметь метод, который мог бы использовать это преимущество.   -  person horchler    schedule 04.09.2013
comment
100x100 крошечный для численных алгоритмов, даже для решения собственных задач. В любом случае, я ничего не знаю о вашей конкретной проблеме. Увеличение ваших дефектных матриц - это не то, чем я когда-либо занимался, поэтому я не могу помочь.   -  person David Heffernan    schedule 04.09.2013
comment
См. мой ответ здесь о том, как символически получить обобщенные собственные векторы для матриц размером более 82 на 82 (предел для моего теста матрица в этом вопросе).   -  person horchler    schedule 14.02.2014
comment
просто прочитал этот пост и проголосовал за повторное открытие вопроса...   -  person bla    schedule 12.02.2015
comment
Спасибо @бла. Дайте мне знать, если у вас есть конкретный интерес к этому. Я мог бы задать вариант этого на SciComp.StackExchange или Math .StackExchange. Я также обновил свой ответ на связанный вопрос — функция jordan по-прежнему ограничена для R2014b. Я не могу комментировать, как это может улучшиться в будущем выпуске... :-)   -  person horchler    schedule 12.02.2015
comment
Просто из интереса, какие матрицы вы пытаетесь здесь диагонализовать? Потому что, если бы они были получены из реальных данных, мне кажется, что они никогда не были бы точно дефектными... Или это должно быть связано с каким-то свойством изучаемой вами системы, и они всегда были бы дефектными точно таким же образом. , поэтому вы должны иметь возможность справляться с этим по-другому... Потому что я не думаю, что есть численно стабильные алгоритмы формы jordan   -  person reverse_engineer    schedule 13.02.2015
comment
@reverse_engineer: Мои матрицы не получены из данных реального слова. Они возникают в результате линеаризации системных ОДУ (конкурентных уравнений Лотки-Вольтерры), параметры которых были специально разработан для того, чтобы якобиан был дефектным (т. е. повторяющимися собственными значениями). Это делается для уменьшения количества свободных параметров и получения желаемого поведения. Пример 3 на 3, приведенный в этом вопросе, является простым случаем. См. мой ответ здесь, чтобы узнать, как создавать более крупные версии этой же системы.   -  person horchler    schedule 13.02.2015
comment
@reverse_engineer: Да, насколько я знаю, для формы Жордана не существует численно устойчивого алгоритма. Но меня не интересует сама форма Жордана, а просто стабильное средство прямого получения обобщенных собственных значений, позволяющее сформировать полный базис — или они неразрывно связаны?   -  person horchler    schedule 13.02.2015
comment
@horchler: Значит, использование чего-то вместе с null((A-eye(3)*eigenvalue)^multiplicity) в соответствии с определением обобщенных собственных векторов вас не удовлетворит?   -  person knedlsepp    schedule 15.02.2015
comment
@horchler: О, неважно! Я только что видел связанный вопрос с этим подходом.   -  person knedlsepp    schedule 15.02.2015
comment
@horchler Да, они связаны, потому что, как только у вас есть обобщенные собственные значения и собственные векторы, ваша форма Жордана определена, или когда у вас есть форма Жордана, вы можете просто прочитать обобщенные собственные векторы ... Theres эта статья, в которой утверждается, что имеется численный алгоритм для вычисления иорданских форм в целом, но насколько я знаю, это нигде не реализовано ... Но если ваши матрицы всегда дефектны точно так же, вы можете обойти это с помощью псевдоинверсий, которые я покажу вам в ответе.   -  person reverse_engineer    schedule 16.02.2015
comment
@reverse_engineer: Я думал об этом. Кажется, правильнее было бы сказать, что алгоритм нахождения обобщенных собственных векторов численно нестабилен? Хотя часть J=V\A*V тоже должна внести свой вклад. Я также видел этот документ/презентацию. Я бы доверял ему немного больше, если бы он был в какой-то степени рецензирован, но это все равно интересно. Я не уверен, что они говорят, что могут решать произвольные системы.   -  person horchler    schedule 17.02.2015


Ответы (1)


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

% example matrix A:
A = [1 0 0 0 0; 
     3 1 0 0 0; 
     6 3 2 0 0; 
     10 6 3 2 0;
     15 10 6 3 2]
% Produce eigenvalues and eigenvectors (not generalized ones)
[vecs,vals] = eig(A)

Это должно вывести:

vecs =

     0         0         0         0    0.0000
     0         0         0    0.2236   -0.2236
     0         0    0.0000   -0.6708    0.6708
     0    0.0000   -0.0000    0.6708   -0.6708
1.0000   -1.0000    1.0000   -0.2236    0.2236


vals =

 2     0     0     0     0
 0     2     0     0     0
 0     0     2     0     0
 0     0     0     1     0
 0     0     0     0     1

Где мы видим, что первые три собственных вектора почти идентичны рабочей точности, как и два последних. Здесь вы должны знать структуру вашей задачи и определить идентичные собственные векторы идентичных собственных значений. Здесь собственные значения абсолютно идентичны, поэтому мы знаем, какие из них рассматривать, и будем считать, что соответствующие векторы 1-2-3 и векторы 4-5 идентичны. (На практике вы, скорее всего, проверите норму разностей собственных векторов и сравните ее с вашим допуском)

Теперь мы приступаем к вычислению обобщенных собственных векторов, но это плохо обусловлено, чтобы решить просто с помощью \ Matlab, потому что, очевидно, (A - lambda*I) не является полным рангом. Итак, мы используем псевдоинверсию:

genvec21 = pinv(A - vals(1,1)*eye(size(A)))*vecs(:,1);
genvec22 = pinv(A - vals(1,1)*eye(size(A)))*genvec21;
genvec1 = pinv(A - vals(4,4)*eye(size(A)))*vecs(:,4);

Что должно дать:

genvec21 =

   -0.0000
    0.0000
   -0.0000
    0.3333
         0

genvec22 =

    0.0000
   -0.0000
    0.1111
   -0.2222
         0

genvec1 =

    0.0745
   -0.8832
    1.5317
    0.6298
   -3.5889

Каковы наши другие обобщенные собственные векторы. Если мы теперь проверим их, чтобы получить нормальную форму Джордана, например:

jordanJ = [vecs(:,1) genvec21 genvec22 vecs(:,4) genvec1];
jordanJ^-1*A*jordanJ

Мы получаем:

ans =

2.0000    1.0000    0.0000   -0.0000   -0.0000
     0    2.0000    1.0000   -0.0000   -0.0000
     0    0.0000    2.0000    0.0000   -0.0000
     0    0.0000    0.0000    1.0000    1.0000
     0    0.0000    0.0000   -0.0000    1.0000

Что является нашей нормальной формой Жордана (с ошибками рабочей точности).

person reverse_engineer    schedule 16.02.2015
comment
Я не уверен, как это может быть численно стабильным? Похоже, что это форма основного метода нахождения обобщенных собственных значений. Для матрицы любого значительного размера со многими повторяющимися собственными значениями можно эффективно умножить pinv(A - vals(1,1)*eye(size(A))) на себя много раз. - person horchler; 17.02.2015
comment
Мне нравится идея каким-то образом каким-то образом воспользоваться структурой - по крайней мере, в некоторых случаях, но это не похоже на то, что вы здесь делаете. Кажется, что можно было бы придумать что-то подобное для моего примера здесь, который имеет только два уникальных собственных значения, одно из которых повторяется N-1 раз. Могут ли в некоторых случаях работать перестановка/сдвиг значений собственного вектора? - person horchler; 17.02.2015
comment
@horchler Эхх, да, возможно, вы правы в том смысле, что если одно собственное значение повторяется большое количество раз, может быть нестабильность этого умножения, но я подозреваю, что на практике это редкость, нет? ... Я думал, вы просто нужен был несимволический способ вычисления обобщенных собственных векторов, и это действительно самый простой способ. Но чем больше я читаю об этом сейчас (вероятно, не так много, как вы читали об этом), тем более очевидным кажется, что этой проблемы избегают и вместо нее используют альтернативы, такие как Schur. В Интернете очень мало информации по этому поводу... В любом случае, удачи в ваших исследованиях! - person reverse_engineer; 17.02.2015