Намиране на обобщени собствени вектори числено в Matlab

Имам матрица като този пример (действителните ми матрици могат да бъдат много по-големи)

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

който има само две линейно независими собствени стойности (собствената стойност -1 се повтаря). Бих искал да получа пълна база с обобщени собствени вектори. Един от начините, по които знам как да направя това, е с функцията jordan на Matlab в 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: Да, видях този въпрос/отговор. Това не е дубликат. Това е основно ръчна процедура, която повтаря това, което е показано в много учебници и в Уикипедия. Знам как сам да получа обобщените собствени стойности - търся нещо малко по-самостоятелно и което използва вградени рутинни процедури, за да внимавам с числените проблеми. Знаете ли нещо или някакви алгоритми? Също така, FYI, 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
Благодаря ти @bla. Кажете ми, ако имате конкретен интерес към това. Може да попитам вариант на това в SciComp.StackExchange или Math .StackExchange. Актуализирах също моя отговор на свързан въпрос – функцията jordan все още е ограничена за R2014b. Все пак не мога да коментирам как може да се подобри в бъдеща версия... :-)   -  person horchler    schedule 12.02.2015
comment
Просто от интерес, какви матрици се опитвате да диагонализирате тук? Защото, ако са извлечени от данни от реалния свят, струва ми се, че никога няма да са точно дефектни... Или трябва да се дължи на някакво свойство на системата, която изучавате, и те винаги ще бъдат дефектни по абсолютно същия начин , така че трябва да можете да се справите с него по различен начин... Защото не мисля, че има числено стабилни алгоритми за йорданска форма   -  person reverse_engineer    schedule 13.02.2015
comment
@reverse_engineer: Моите матрици не са извлечени от реални думи. Те идват от линеаризиране на система ODE (конкурентни уравнения на Лотка-Волтера), чиито параметри са били специално проектирани така, че Якобианът да е дефектен (т.е. повтарящи се собствени стойности). Това се прави, за да се намали броят на свободните параметри и да се получи желано поведение. Примерът 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 Да, двете са свързани, защото след като имате обобщени собствени стойности и собствени вектори, вашата Йорданова форма е дефинирана, или когато имате Йорданова форма, можете просто да прочетете обобщените собствени вектори... Има тази статия, която твърди, че разполага с числено алго за изчисляване на jordan форми като цяло, но не е приложен никъде, доколкото знам... Но ако вашите матрици винаги са дефектни по абсолютно същия начин, можете да го заобиколите с псевдообратни, които ще ви покажа в отговор.   -  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