Vettori sparsi

[Studente Anonimo]

Come faccio a creare dei vettori sparsi su MatLab?
Dove la prima componente e l'ultima sono uguali a 1 e tutte le altre componenti sono uguali a 0.
Il problema con il mio codice è che ho matrici e vettori che vanno fino a dimensione \( 10^7 \times 10^7 \) e \( 1 \times 10^7 \), la matrice sono riuscito a crearla sparsa i vettori no, quindi è un casino! Perché quando la dimensione arriva a \(10^5 \) mi occupa 74 GB di memoria e mi blocca l'algoritmo....

[feddy]

Non uso MatLab da un po', ma siccome lì "tutto è matrice", potresti provare a definire ad utilizzare sparse anche sul vettore, definendolo come una matrice con una riga (o una colonna).

Qual è l'obiettivo del tuo codice? Con un po' di contesto probabilmente è possibile esserti più utile

[Studente Anonimo]

Devo creare un algoritmo che richieda un numero al più \(O(n)\) di operazioni per risolvere un sistema
\( Bx = b \)
dove \( b \) è un vettore solo fatto di \(1\) e \(B\) è una matrice della forma
\[ B = \begin{pmatrix}
-1& 1 & & & 1 \\
1& -2 & \ddots & & \\
& \ddots & \ddots & \ddots & \\
& & \ddots & -2 & 1\\
1& & & 1 & -1
\end{pmatrix} \]
di dimension \(n \times n \) con \( n=10,10^2,\ldots,10^7 \).
Quello che ho fatto io è stato creare una matrice tridiagonale
\[ A = \begin{pmatrix}
-2& 1 & & & 0 \\
1& -2 & \ddots & & \\
& \ddots & \ddots & \ddots & \\
& & \ddots & -2 & 1\\
0& & & 1 & -2
\end{pmatrix} \]
e due vettori \( u=v=(1,0,\ldots,0,1) \) in modo tale che posso fare
\[ B = A + uv^{T} \]
e grazie a Sherman-Morrison so l'inversa di \( B \) facilmente inoltre siccome \(A\) è tridiagonale abbiamo un algoritmo che è \(O(n) \). Per risolvere \( Ax = b \).

[Studente Anonimo]

L'idea del mio codice è questa:

function alg = my_alg(n)
b = ones(10^n,1);
A=sparse_matrix(10^n);
u = zeros(1,10^n);
v = zeros(10^n,1);
u(10^n)=1;
u(1)=1;
%u;
v(10^n)=1;
v(1)=1;
sparse(v);
sparse(u);
C = v*u;
B= A + C;
x1 = A \ b ;
x2 = A \ v;
x = x1 - ( x2*u*x1 )./(1+u*x2);
alg = x;
end

ho provato anche a fare direttamente

function alg = my_alg(n)
b = ones(10^n,1);
A=sparse_matrix(10^n);
u = zeros(1,10^n);
v = zeros(10^n,1);
u(10^n)=1;
u(1)=1;
%u;
v(10^n)=1;
v(1)=1;
sparse(v);
sparse(u);
C = v*u;
B= A + C;
x= B\ b;
alg = x;
end

ma in entrambi i casi mi dice la stessa cosa: che ho superato la dimensione per array.
ps: ho creato una funzione così posso dopo verificare con timeit che sia effettivamente un tempo O(n).

[Studente Anonimo]

Digitando "help sparse" sul terminale mi dice questo: ma non capisco come cavolo fare i due vettori di cui ho bisogno, perché non capisco come dichiararli onestamente.

sparse Create sparse matrix.
S = sparse(X) converts a sparse or full matrix to sparse form by squeezing
out any zero elements.

S = sparse(i,j,s,m,n,nzmax) uses vectors i, j, and s to generate an
m-by-n sparse matrix such that S(i(k),j(k)) = s(k), with space
allocated for nzmax nonzeros. Vectors i, j, and s are all the same
length. Any elements of s that are zero are ignored, along with the
corresponding values of i and j. Any elements of s that have duplicate
values of i and j are added together. The argument s and one of the
arguments i or j may be scalars, in which case the scalars are expanded
so that the first three arguments all have the same length.

S = sparse(i,j,s,m,n) where nzmax = length(s).

S = sparse(i,j,s) where m = max(i) and n = max(j).

S = sparse(m,n) abbreviates sparse([],[],[],m,n,0). This
generates the ultimate sparse matrix, an m-by-n all zero matrix.

For example, this dissects and then reassembles a sparse matrix:

[i,j,s] = find(S);
[m,n] = size(S);
S = sparse(i,j,s,m,n);

So does this, if the last row and column have nonzero entries:

[i,j,s] = find(S);
S = sparse(i,j,s);

[Studente Anonimo]

So.. scemo però... non dichiaravo una nuova variabile a sparse(v) e sparse(u) e non è che mi trasforma v e u in sparse... e quindi quando usavo v e u sotto non me le considerava sparse.

[feddy]

Ho visto ora i messaggi. Capita a tutti di fare errori così

Per quanto riguarda l'esercizio in sè... Sherman-Morrison è utile, ma ci sono un po' di sistemi lineari (chiaramente l'inversa da sola non si costruisce) da risolvere usando quella formula, e quelli sono il collo di bottiglia in termini di complessità, ovviamente. Non hai detto come li risolvi

[feddy]

Ok ho visto il codice, usi backslash. Per sistemi sparsi, di grandi dimensioni, il backslash di MatLab può essere altamente inefficiente. Tuttavia nel tuo caso funziona ancora tutto quanto bene.

Un altro modo, più standard, per risolvere in $O(n)$ un sistema di quel tipo è con l'algoritmo di Thomas. In pratica altro non è che eliminazione di Gauss applicata nel caso in cui il sistema sia tridiagonale (il fatto che ci siano quei due $1$ in alto e in basso non complica le cose più di tanto). E' un algoritmo standard e su google trovi già tutto quello che ti serve. Nel caso tu lo voglia implementare, puoi notare come non ti serve storare alcuna matrice, ma solamente le entrate non nulle in vettori e poi operare su di essi.

Quella matrice tra l'altro è esattamente la discretizzazione con differenze finite di $u''$ con condizione periodiche (a meno di un fattore $h^{-2}$).

[axpgn]

@feddy
[ot]Tra "storare" e "mecciano" sei proprio un informatico doc [/ot]


Cordialmente, Alex

[feddy]

@Alex

[ot]AHAHA usando quei termini in effetti passo proprio per un informatico, ma in realtà sono un numerico

Hai ragione però, in effetti storare non è italiano... meglio usare "allocare" [/ot]

[axpgn]

@feddy
[ot]"storare" già era standard quarant'anni fa , e "mecciano" giù di lì pero quest'ultimo è proprio da informatici "puri" mentre l'altro ha un'audience più ampia (ci ho messo un po' a capire il significato del secondo mentre "storà" dalle mie parti significa "stancare, stufare" ma più significativamente )[/ot]

[feddy]

[ot]Oh cielo, non sapevo che "mecciano" si dicesse, io l'avevo scritto per scherzare[/ot]

[axpgn]

Oh, sì

Anzi, si declina pure