Despre descopunerea in vectori proprii si descompunerea SVD

Sa tin minte, ca m-oi face o data de rusine prin vreun articol:

SVD:
- valori singulare reale si pozitive
- U si V ortogonale, dar distincte

Eig:
- valori proprii negative sau chiar si complexe
- V poate for oricum, nu neaparat ortogonala (dar e una singura, V si V^-1)

La matrici simetrice:
-  vectorii proprii sunt ortogonali, V = ortogonala, dar in continuare e una singura, V si V', iar valorile proprii pot fi negative

Asadar, in cazul matricilor reale simetrice, cele doua NU SUNT echivalente (ar trebui sa fie si SPD pentru asta, ceea ce oricum nu poti afla decat facand eig):
- M = U D U'
- M = U S V'

DAR coloanele din U si V pot avea semnul diferit (in schimb S sunt musai pozitive), in timp ce coloanele lui U si U^-1=U' au acelasi semn (in schimb semsnul se muta la valorile proprii, care pot fi negative).

Asadar poti sa gasesti SVD din EIG in acest caz, daca muti semnul minus de la valorile proprii negative (daca sunt) la unul din vectorii din U sau V, stricand astfel simetria dintre U si U', care devin U si V. Trebuie asadar sa folosesti eigs(largest magnitude), apoi muti minusul valorilor proprii negative la V.

Nearest low-rank, apoi extras square root dintr-o matrice simetrica:

Low-rank faci cu SVD, care muta semnul minus de la valorile proprii la V, stricand simetria dintre U si V. Dar ca sa extragi apoi square root, adica M = X'*X, ai nevoie de U si V simetrice si valori proprii nenule, simultan (pentru cu X=USV' implica  X'X = VS'SV = VS^2V' = semipozitiv definita). Asadar tre sa muti minusurile invers de la V la valorile proprii, apoi retii doar valorile proprii pozitive. Deci tot ramai doar cu valorile proprii largest amplitude, pozitive. Sau cel mai bine, luand ambele etape in considerare, faci de la bun inceput eigs(largest amplitude) si pastrezi cele mai mari pozitive.

Produsul a doua matrici pozitiv definite este pozitiv definit?

Raspuns: Nu neaparat

Interpretare geometrica:

Pozitiv definit inseamna:
$x^T A x >0, \forall x,$
care poate fi rescris ca produsul scalar intre $x$ si $Ax$
$x^T A x  = \langle Ax, x\rangle >0$
Asta inseamna ca vectorul imagine $Ax$ este in acelasi semiplan ca si vectorul $x$ (la dreapta liniei segmentate in figura, unde vectorul $x = OA$).

Acum, matricea B*A este pozitiv definita daca
$x^T BA x > 0, \forall x,$
care poate fi rescris ca produsul scalar intre $Ax$ si $B^Tx$:
$x^T BA x  = \langle B^Tx, Ax \rangle > 0, \forall x.$

In figura avem insa un vectorii $Ax$ (OC) si $B^Tx$ (OD, din greseala este notat Bx) care sunt amandoi in acelasi semiplan cu $x$, dar ei intre ei formeaza un unghi mai mare de 90 grade si deci produsul scalar este negativ. Deci $x^T BA x < 0$.

Importanta simetriei:
O matrice pozitiv definita are toate valorile proprii mai mari ca zero. Dar reciproca nu e intotdeauna valabila, cum se vede in cele de mai jos. O matrice nesimetrica poate avea valori proprii pozitive dar totusi sa nu fie pozitiv definita. Doar daca mai este si simetrica atunci se poate concluziona ca matrice simetrica + valori proprii > 0 implica pozitiv defninita.
A se vedea referinta: http://math.stackexchange.com/questions/4336/if-eigenvalues-are-positive-is-the-matrix-positive-definite


Concluzie: testarea valorilor proprii nu este suficienta pentru a garanta pozitiv definitia decat daca matricea este si simetrica.

Script Matlab:
% Doua matrici pos.def. al caror produs nu e pos.def.:

A = diag([0.1, 100]);
B = [1 -1; 1 1]^-1 * diag([10 0.1]) * [1 -1; 1 1];
x = [1; 0.1];

% Simetrice si eigs > 0 <=> pos.def.
eig(A)
eig(B)

% B*A nu e posdef pentru ca x'(BA)x < 0
x' * B*A * x

% desi are eigs > 0
eig(B*A)

% pentru ca lipseste conditia de simetrie
% Vezi si
% http://math.stackexchange.com/questions/4336/if-eigenvalues-are-positive-is-the-matrix-positive-definite

Dezactivare system beep sub Linux (non-root)

xset -b b 0

Extrem de folositor atunci cand Matlab-ul te innebuneste cu beep-uri in editor si la erori.

Aventuri cu Matlab mex sub linux x86_64

Am un script Matab care foloseste bucata de cod de pe net care foloseste un third-party mexeig.c, care la randul sau foloseste functia dsyevd() din Lapack. Sub Windows, mai demult nu mai stiu cum, compilasem .c -ul in mex si mergea foarte bine.Din pacate, intreg script-ul dura vreo doua zile pe 3 threaduri ale PC-ului de la birou.

Azi, ca sa nu imi ocup din nou comp-ul pentru 2 zile, mi-am zis sa copii codul pe serverul sub Linux si sa-l rulez acolo pe 6 threaudri (tip: startworkers 6, si apoi matlabpool serverplus).

Cand colo, pac eroare ca nu stie mexeig(). Dupa ca ma gandesc putin, imi dau seama ca nu compilasem mexeig.c sub noua arhitectura. Bine, bag repede un mex mexeig.c. Pac, nu gaseste dsyevd la link-editare, undefined reference. Bine din nou, bag reloaded un mex mexeig.c -llapack. Pac, skipping incompatible library vreo 4 lapack-uri si in continuare undefined reference. Si in plus, la deruta, si un warning ca folosesc gcc 3.4.6 in timp ce minimul testat de ei este 4.0.0.

Eu ma gandesc, bineinteles, ca e problema cu gcc-ul. Injur in gand administratorul care tine o versiune antica si acceptata de CentOS 4.8, cu un gcc antediluvian, si ma gandesc sa-mi instalez o varianta locala de gcc 4.x. Pierd vreo 40 de minute pe net, documentatie, cauta repository CentOS 5.x, trag un rpm cu gcc-4.1.2, ii zic de la obraz un rpm -i --test gcc..., bineintelesc ca nu are dependecies necesare. Mama lui de administrator, imi zic in gand. Dead end, timp pierdut.

Apoi, ca din intamplare, bag un gcc4. De fapt nu din intamplare, vazusem prin lib un folder cu niste librarii gcc 4.1.2. Ce sa vezi, gcc4 merge. E instalat si gcc 4.1.2, dar sub numele de gcc4, nu gcc. Bine, zic, acum cum dracu' il fac sa-l vada si matlab-ul, ca in mex -setup imi apare doar gcc-ul vechi. Iar documentatie, aflu ca exista un mexopts.sh care defineste comenzile si optiunile folosite de Matlab la compilare, intru, replace gcc cu gcc4, dau victorios un mex mexeig.c -llapack, la fel. Skipping incompatible library, undefined reference. Mama naibii. Macar warning--ul se dusese, deci folosea versiune noua de gcc. Toate astea au durat inca vreo juma de ora.

Exasperat, dau un View (F3) din Midnight Commander pe liblapack.so.3.0.3, sa vad, ce dracu' nu gaseste dsyevd, te pomeni ca nu e definita in librarie, sau mai stiu io ce. Si-mi sare in ochi un i386. Nu realizez pe moment, dar dupa vreo 5-10 minute subconstientul ma trage de maneca: i386 ?!? Pe 64 biti?

Mama ei de porcarie, exista un folder lib64 in usr, dar acolo n-are lapack. Asta e, imi zic, CentOS-ul asta are numai versiunea perfect inutila i386, de ce naiba n-are si pe 64 de biti? Iar fugi pe net, cauta un repository, trag un rpm cu liblapack x86_64, dezarhiveaza intr-un folder local, mex -llapack -L/home/... mexeig.c. L-am prins, zic. Cand colo, stupoare. Undefined reference. Nu mai da skipping incompatible library, deci libraria e buna, de ce naiba nu merge?

Ma uit cu F3 in librarie, caut dsyevd, exista, dar vad si multe _ si imi aduc aminte ca e o chestie cu astea. Intru in mexeig.c si inlocuiesc apelul lui dsyevd cu dsyevd_ (imi aduc aminte ca am mai facut asta odata, probabil cand am compilat pe Windows). In sfarsit, merge!

In total, am irosit vreo ora si jumatate, plus un sfert de ora postul asta.

Lectii din treaba asta:
1. Skipping incompatible library se refera la arhitectura pe care a fost compilata libraria (i386, x86_64 etc). nu la versiunea gcc sau mai stiu eu ce.
2. Poti sa adaugi un _ la numele functiei cand nu merge sa link-editezi cu librarii din astea publice si n-ai idee de ce.
3. gcc mai poate avea si numele gcc4

Update. Ehei, dragii mosului, eu inca nu testasem daca mexeig-ul compilat chiar merge. Adica m-am oprit la compilarea si link-area cu succes. Mai departe, cand dau sa rulez, sare cu o eroare run-time cum ca nu gaseste liblapack-ul. Noroc ca vazusem in help-ul lui mex ca, daca compilez cu o librarie dintr-un folder oarecare, tre sa modific si path-urile ca Matlab-ul chiar sa gaseasca libraria respectiva la run-time. Asa ca le-am adaugat frumusel intr-un script fix inainte de a starta matlab-ul (matlab-ul fiind pornit din script, e child process si mosteneste environmentul acestuia):

#!/bin/bash

LD_LIBRARY_PATH=:/home/...$LD_LIBRARY_PATH
export LD_LIBRARY_PATH

cd "CS/Ale mele"
matlab

Au trebuit cea modificari fata de ce zicea help-ul Matlabului, adica am pus ":" inainte de calea mea, si nu intre (pentru ca LD_L... default continea un ":" la inceput), si a trebuit si "$", desigur.

Si inca nu m-a lasat. Cica nu gaseste nici libraria blas. Asa ca ia-o si pe asta de pe net si dezarhiveaz-o in acelasi usr local. Intr-un final, merge.

Update 2  Merge pe dracu! Da, merge daca rulez pe un singur thread. Daca insa incep un matlabpool si rulez codul in parfor, s-a dus. Iarasi nu gaseste liblapack-ul. Problema este ca, in acest caz, codul este executat de catre workers, care nu mostenesc environmentul de la Matlab-ul mare, pentru ca ei sunt creati independent, cu daemon-ul mdce. Asa ca ce-am facut: in primul rand, am adaugat noua LD_LIBRARY_PATH atat in bashrc cat si in bash_profile, si le-am scos din script-ul meu. Sa fie mereu definit mereu si peste tot, ce naiba. Apoi am dat un source, cica, apoi am si repornit vncserver ca sa ma loghez inca o data (nu stiu daca chiar a fost un logout pe bune) ca sa fiu sigur ca s-au citit fisierele din nou. Apoi am reastartat daemon-ul mdce. Asta era chestia esentiala. Acum mdce a re-pornit dintr-un terminal, care mostenise environment-ul bun de la bashrc, deci si mdce il mosteneste, si apoi toti workerii creati de mdce. Se pare ca, intr-un sfarsit, merge (a treia oara)!