Generuj symetryczną dodatnią określoną macierz z wcześniej określonym wzorcem sparsity

9

Próbuję wygenerować macierz korelacji p×p (symetryczny psd) z wcześniej określoną strukturą sparsity (określoną przez wykres na pwęzły). Węzły połączone na wykresie mają korelacjęρU(0,1), reszta to 0, a przekątna to 1.

Próbowałem wygenerować tę macierz kilka razy, ale rzadko otrzymuję prawidłową macierz korelacji.

Czy istnieje sposób, aby zapewnić whp macierzy korelacji? Zauważ, że mogę mieć tylko dodatnią korelacjęρU(1,1) itp. nie jest opcją.

Każda pomoc jest mile widziana!

Łowca ostrzy
źródło
Może funkcja nearPD pakietu Matrix w R może pomóc.
niandra82
Jaka jest twoja miara rzadkości, która jest dla ciebie ustalona? Czy twoje dane powinny być ciągłe binarne czy nieujemne?
ttnphns
@ niandra82: nearPD nie jest dobre, ponieważ zniszczy rzadkość matrycy.
Blade Runner
1
Zasadniczo nie ma takich rozkładów macierzy, jak opisano w tym pytaniu. Rozważmy na przykład3×3 skrzynia z trzema współczynnikami ρ,σ,τ. Gdybyτ=0 i ρ>0,σ>0, następnie ρ2+σ2<1wtedy i tylko wtedy, gdy macierz jest pozytywnie określona. Ale wtedy nie możesz mieć obuρU(0,1) i σU(0,1).
whuber
3
To dlaczego nie najpierw wygenerować macierz korelacji. Następnie utwórz indeks symetryczny dla tej macierzy, w której zmusisz indeksowane elementy do 0. Niewielkość byłaby określona przez rozmiar indeksu, a możesz wstawić losowość za pomocą funkcji podobnej do próbki w r. Bez względu na to, ile elementów przekątnych wymusisz na 0, matix nadal będzie pd
Zachary Blumenfeld

Odpowiedzi:

2

Zamknij, ale nie cygaro dla @Rodrigo de Azevedo.

Rozwiązaniem jest użycie programowania półfinałowego, aby znaleźć maksymalną wartość, ρmaxoraz minimalna wartość (z zastrzeżeniem nieujemności), ρmin, z ρtak, że macierz korelacji z zalecanym wzorem rzadkości jest dodatnia półfinałowa (psd). Wszystkie wartościρ takie, że ρmaxρρmax, stworzy macierze psd (ćwiczenie dla czytelnika)

Dlatego musisz albo wybrać rozkład ρ który może przyjmować tylko wartości [ρmax,ρmax], lub musisz użyć akceptacji / odrzucenia i odrzucić wszelkie wygenerowane wartości ρ które nie produkują macierzy psd.

Przykład matrycy 4 na 4 przy użyciu YALMIP w MATLAB

sdpvar rho % declare rho to be a scalar variable
% find maximum value of rho (by minimizing -rho) subject to prescribed matrix being psd.
optimize([1 0 rho 0;0 1 rho 0;rho rho 1 rho;0 0 rho 1] >= 0,-rho) 
% find minimum value of rho subject to prescribed matrix being psd and rho being >= 0.
optimize([[1 0 rho 0;0 1 rho 0;rho rho 1 rho;0 0 rho 1] >= 0,rho >= 0],rho) 

Wyniki: maksymalne rho = 0,57735, minimalne rho = 0. Zrozumiałe jest, że zero będzie minimalną wartością rho pod warunkiem, że rho będzie nieujemne, a zalecaną macierzą będzie psd, niezależnie od wzorca wymiaru lub rzadkości. Dlatego nie jest konieczne uruchamianie optymalizacji półfinałowej, aby znaleźć minimalną nieujemną wartośćρ.

Mark L. Stone
źródło
4
Jest to interesująca interpretacja pytania: zakłada, że ​​wszystkie niezerowe współczynniki off-diagonalne są równe (tym samym ogromnie upraszczając problem). Nie jest jasne, czy taka była zamierzona interpretacja, czy też wszystkie niezerowe współczynniki off-diagonalne mają być niezależnymi realizacjami od wspólnego rozkładu.
whuber
Taką interpretację zrobiłem. Teraz, kiedy o tym wspominasz, mogłem zobaczyć inną interpretację. Przynajmniej moja interpretacja ma tę zaletę, że powoduje dość dobrze zdefiniowany problem. Przypuszczam, że można sformułować problem, którego rozwiązania nie zbadałem, aby znaleźć maksymalną wartość ρ w taki sposób, że wszystkie niezerowe elementy o przekątnej jednego trójkąta macierzy korelacji mogą być wypełnione niekoniecznie równymi nieujemnymi wartościami ≤ tę wartość i koniecznie sprawi, że w pełni wypełniona macierz będzie psd.
Mark L. Stone,
0

Macierz korelacji jest symetryczna, dodatnia półfinałowa i ma 1jest na swojej głównej przekątnej. Można znaleźćn×nmacierz korelacji poprzez rozwiązanie następującego programu półfinałowego (SDP), w którym funkcja celu jest dowolna, powiedzmy, funkcja zerowa

minimizeOn,Xsubject tox11=x22==xnn=1XOn

Jeśli ktoś ma dodatkowe ograniczenia, takie jak ograniczenia rzadkości

xij=0 for all (i,j)Z[n]×[n]

i ograniczenia nieujemności, XOn, a następnie rozwiązuje następujący SDP

minimizeOn,Xsubject tox11=x22==xnn=1xij=0 for all (i,j)Z[n]×[n]XOnXOn

ZA 3×3 przykład

Załóżmy, że chcemy mieć x13=0 i x12,x230. Oto skrypt MATLAB + CVX ,

cvx_begin sdp

    variable X(3,3) symmetric

    minimize( trace(zeros(3,3)*X) )
    subject to

        % put ones on the main diagonal
        X(1,1)==1
        X(2,2)==1
        X(3,3)==1

        % put a zero in the northeast and southwest corners
        X(1,3)==0

        % impose nonnegativity
        X(1,2)>=0
        X(2,3)>=0

        % impose positive semidefiniteness
        X >= 0

cvx_end

Uruchamianie skryptu

Calling sedumi: 8 variables, 6 equality constraints
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 6, order n = 6, dim = 12, blocks = 2
nnz(A) = 8 + 0, nnz(ADA) = 36, nnz(L) = 21
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            3.00E+000 0.000
  1 : -1.18E-001 6.45E-001 0.000 0.2150 0.9000 0.9000   1.86  1  1  1.2E+000
  2 : -6.89E-004 2.25E-002 0.000 0.0349 0.9900 0.9900   1.52  1  1  3.5E-001
  3 : -6.48E-009 9.72E-007 0.097 0.0000 1.0000 1.0000   1.01  1  1  3.8E-006
  4 : -3.05E-010 2.15E-009 0.000 0.0022 0.9990 0.9990   1.00  1  1  1.5E-007
  5 : -2.93E-016 5.06E-015 0.000 0.0000 1.0000 1.0000   1.00  1  1  3.2E-013

iter seconds digits       c*x               b*y
  5      0.3   5.8  0.0000000000e+000 -2.9302886987e-016
|Ax-b| =  1.7e-015, [Ay-c]_+ =  6.1E-016, |x|= 2.0e+000, |y|= 1.5e-015

Detailed timing (sec)
   Pre          IPM          Post
1.563E-001    2.500E-001    1.094E-001    
Max-norms: ||b||=1, ||c|| = 0,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 1.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +0

Zobaczmy, jakie rozwiązanie znalazło CVX,

>> X

X =

    1.0000    0.4143         0
    0.4143    1.0000    0.4143
         0    0.4143    1.0000

Czy ta matryca jest pozytywna? Pozytywne określone?

>> rank(X)

ans =

     3

>> eigs(X)

ans =

    1.5860
    1.0000
    0.4140

Jest pozytywny określony, zgodnie z oczekiwaniami. Korzystając z niezerowej (liniowej) funkcji celu, możemy znaleźć dodatnie macierze korelacji półfinałowej.

Rodrigo de Azevedo
źródło
Ponieważ na tej stronie „generowanie” będzie rozumiane jako „losowanie z losowego rozkładu”, czy mógłbyś wyjaśnić, w jaki sposób twój kod generuje macierze losowej korelacji i wskazać, jaki rozkład następuje?
whuber
@whuber OP prosi o niemożliwe. Skomentowałeś to 1 stycznia 2015 r. Jeśli chcesz wygenerować losowe macierze korelacji, to wygeneruj losową macierz kwadratową i użyj jej w funkcji celu w powyższym programie półfinałowym. Lub wygeneruj realizacje zmiennej losowej, która jest równomierna na kostce umieść je w wpisach macierzy (korelacji) z na główną przekątną i odrzuć te, które nie są dodatnie w połowie skończone. Jeśli istnieją ograniczenia nieujemności, równomiernie próbkuj kostkę
[1,1](n2)
1
[0,1](n2)
Rodrigo de Azevedo
3
@whuber Oto trójwymiarowy eliptyczny [png], który odwzorowuje na zestaw macierzy korelacji . To, czego chce OP, to przecięcie eliptonu z nieujemną oktantą, a następnie przecięcie go z płaszczyznami postaci . Jeśli macierz ma , to musi znajdować się we wnętrzu eliptopu. Korzystając z SDP z niezerowymi funkcjami celu, można próbkować powierzchnię eliptyczna. Ponieważ eliptop jest wypukły, wypukłe kombinacje punktów powierzchni będą również mapowane na macierze korelacji. 3×3xij=00
Rodrigo de Azevedo
1
To doskonały sposób na opisanie sytuacji.
whuber
3
Masz rację co do tego, jak zmniejszają się względne objętości. Właśnie dlatego jest to trudny problem.
whuber