Algebra lineare

La traduzione di questa pagina non è aggiornata. Fai clic qui per vedere l'ultima versione in inglese.

Le matrici nell'ambiente di MATLAB

Questo argomento contiene un'introduzione alla creazione di matrici e all'esecuzione di calcoli di base sulle matrici in MATLAB^®.

L'ambiente di MATLAB utilizza il termine matrice per indicare una variabile con numeri reali o complessi disposti in una griglia bidimensionale. Un array è, più in generale, un vettore, una matrice o una griglia di numeri di maggiori dimensioni. Tutti gli array in MATLAB sono rettangolari, nel senso che i vettori che lo compongono, lungo ogni dimensione, presentano tutti la stessa lunghezza. Le operazioni matematiche definite sulle matrici sono oggetto dell'algebra lineare.

Creazione di matrici

MATLAB presenta numerose funzioni che creano tipi diversi di matrici. Ad esempio, è possibile creare una matrice simmetrica con elementi basati sul triangolo di Pascal:

A = pascal(3)

A =
       1     1     1
       1     2     3
       1     3     6

È anche possibile creare una matrice di un quadrato magico asimmetrica, con somme uguali per le righe e le colonne:

B = magic(3)

B =
       8     1     6
       3     5     7
       4     9     2

Un altro esempio è una matrice rettangolare 3x2 di numeri interi casuali. In questo caso il primo input in randi descrive il range di valori possibili per i numeri interi, mentre i due input successivi descrivono il numero delle righe e delle colonne.

C = randi(10,3,2)

Un vettore colonna è una matrice mx1, un vettore riga è una matrice 1xn e uno scalare è una matrice 1x1. Per definire manualmente una matrice, utilizzare le parentesi quadre [ ] per indicare l'inizio e la fine di un array. All'interno delle parentesi utilizzare un punto e virgola ; per indicare la fine di una riga. Gli scalari (matrice 1x1) non richiedono parentesi. Ad esempio, queste dichiarazioni producono un vettore colonna, un vettore riga e uno scalare:

u = [3; 1; 4]

v = [2 0 -1]

s = 7

u =
       3
       1
       4

v =
       2     0    -1

s =
       7

Per maggiori informazioni su come creare matrici e sulle relative operazioni, vedere Creating, Concatenating, and Expanding Matrices.

Somma e sottrazione di matrici

La somma e la sottrazione di matrici e array viene eseguita un elemento per volta, o in modalità orientata verso gli elementi. Ad esempio, aggiungendo A a B e quindi sottraendo A dal risultato, si ottiene B:

X = A + B

X =
       9     2     7
       4     7    10
       5    12     8

Y = X - A

Y =
       8     1     6
       3     5     7
       4     9     2

La somma e la sottrazione richiedono che entrambe le matrici presentino dimensioni compatibili. Se le dimensioni non sono compatibili viene generato un errore:

X = A + C

Error using  + 
Matrix dimensions must agree.

Per maggiori informazioni, vedere Array vs. Matrix Operations.

Prodotti e trasposizioni di vettori

Un vettore riga e un vettore colonna della stessa lunghezza possono essere moltiplicati in qualsiasi ordine. Il risultato è o uno scalare, denominato prodotto interno, o una matrice, denominata prodotto esterno:

u = [3; 1; 4];
v = [2 0 -1];
x = v*u

x =

     2

X = u*v

X =

     6     0    -3
     2     0    -1
     8     0    -4

Per le matrici reali, l'operazione di trasposizione scambia a_ij e a_ji. Per le matrici complesse occorre considerare anche se utilizzare il coniugato complesso degli elementi complessi nell'array per formare il trasposto complesso coniugato. MATLAB utilizza l'operatore apostrofo (') per eseguire un trasposto complesso coniugato e l'operatore punto-apostrofo (.') per il trasporto senza coniugazione. I due operatori restituiscono lo stesso risultato in caso di matrici che contengono tutti elementi reali.

La matrice dell'esempio A = pascal(3) è simmetrica, quindi A' è uguale ad A. B = magic(3) invece non è simmetrica, quindi B' presenta gli elementi riflessi lungo la diagonale principale:

B = magic(3)

B =

     8     1     6
     3     5     7
     4     9     2

X = B'

X =

     8     3     4
     1     5     9
     6     7     2

Per i vettori, la trasposizione trasforma un vettore riga in vettore colonna e viceversa:

Se sia x che y sono vettori colonna reali, il prodotto x*y non è definito, ma i due prodotti

x'*y

y'*x

generano lo stesso risultato scalare. Questa quantità viene utilizzata tanto di frequente da presentare tre nomi diversi: prodotto interno, prodotto scalare o prodotto punto. È disponibile persino una funzione dedicata per i prodotti punto, dal nome dot.

Per le matrici o i vettori complessi z, la quantità z' non traspone solo il vettore o la matrice, ma converte anche ogni elemento complesso nel corrispettivo coniugato complesso. Ciò significa che cambia il segno della parte immaginaria di ogni elemento complesso. Si consideri ad esempio la matrice complessa

z = [1+2i 7-3i 3+4i; 6-2i 9i 4+7i]

z =

   1.0000 + 2.0000i   7.0000 - 3.0000i   3.0000 + 4.0000i
   6.0000 - 2.0000i   0.0000 + 9.0000i   4.0000 + 7.0000i

Il trasposto complesso coniugato di z è:

z'

ans =

   1.0000 - 2.0000i   6.0000 + 2.0000i
   7.0000 + 3.0000i   0.0000 - 9.0000i
   3.0000 - 4.0000i   4.0000 - 7.0000i

Il trasposto complesso non coniugato, in cui viene mantenuto il segno della parte complessa di ciascun elemento, è denotato con z.':

z.'

ans =

   1.0000 + 2.0000i   6.0000 - 2.0000i
   7.0000 - 3.0000i   0.0000 + 9.0000i
   3.0000 + 4.0000i   4.0000 + 7.0000i

Per i vettori complessi, i due prodotti scalari x'*y e y'*x sono coniugati complessi l'uno dell'altro e il prodotto scalare x'*x di un vettore complesso rispetto a se stesso è reale.

Moltiplicazione di matrici

La moltiplicazione di matrici è definita in modo da riflettere la composizione delle trasformazioni lineari sottostanti e da consentire la rappresentazione compatta di sistemi di equazioni lineari simultanee. Il prodotto di matrice C = AB è definito quando la dimensione della colonna di A è uguale alla dimensione della riga di B, oppure quando uno di questi elementi è uno scalare. Se A è mxp e B è pxn, il loro prodotto C è mxn. Il prodotto in realtà può essere definito con i loop MATLAB for, con la notazione colon e con i prodotti punto del vettore:

A = pascal(3);
B = magic(3);
m = 3; 
n = 3;
for i = 1:m
     for j = 1:n
        C(i,j) = A(i,:)*B(:,j);
     end
end

MATLAB utilizza un asterisco per denotare la moltiplicazione di matrici, come in C = A*B. La moltiplicazione di matrici non è commutativa; A*B infatti non è uguale a B*A:

X = A*B

X =
      15    15    15
      26    38    26
      41    70    39

Y = B*A

Y =
      15    28    47
      15    34    60
      15    28    43

Una matrice può essere moltiplicata sulla destra da un vettore colonna e sulla sinistra da un vettore riga:

u = [3; 1; 4];
x = A*u

v = [2 0 -1];
y = v*B

y =

    12    -7    10

Le moltiplicazioni di matrici rettangolari devono soddisfare la condizione di compatibilità delle dimensioni. Poiché A è 3x3 e C è 3x2, è possibile moltiplicarli per ottenere un risultato 3x2, (la dimensione interna comune viene annullata):

X = A*C

La moltiplicazione tuttavia non funziona nell'ordine inverso:

Y = C*A

Error using  * 
Incorrect dimensions for matrix multiplication. Check that the number of columns 
in the first matrix matches the number of rows in the second matrix. To perform 
elementwise multiplication, use '.*'.

È possibile moltiplicare qualsiasi elemento con uno scalare:

s = 10;
w = s*y

w =

   120   -70   100

Quando si moltiplica un array per uno scalare, quest'ultimo si espande implicitamente per presentare le stesse dimensioni dell'altro input. Questo spesso prende il nome di espansione scalare.

Matrici di identità

La notazione matematica generalmente accettata utilizza la lettera maiuscola I per la denotazione di matrici di identità, cioè matrici di varie dimensioni con diversi elementi uno sulla diagonale principale ed elementi zero altrove. La proprietà di queste matrici è: AI = A e IA = A, quando le dimensioni sono compatibili.

La versione originale di MATLAB non consentiva l'uso di I a questo scopo perché non distingueva tra maiuscole e minuscole e i serviva già come pedice e come unità complessa. Il problema è stato risolto con un gioco di parole in lingua inglese fondato sulla medesima pronuncia dei termini I e eye. La funzione

eye(m,n)

restituisce una matrice identità rettangolaremxn, mentre eye(n) restituisce una matrice identità quadratanxn.

Inverso delle matrici

Se una matrice A è quadrata e non singolare (con determinante non zero), le equazioni AX = I e XA = I generano la stessa soluzione X. Questa soluzione è detta inverso di A ed è denotata con A^-1. Sia la funzione inv che l'espressione A^-1 calcolano l'inverso della matrice.

A = pascal(3)

A =
       1     1     1
       1     2     3
       1     3     6

X = inv(A)

X =

    3.0000   -3.0000    1.0000
   -3.0000    5.0000   -2.0000
    1.0000   -2.0000    1.0000

A*X

ans =

    1.0000         0         0
    0.0000    1.0000   -0.0000
   -0.0000    0.0000    1.0000

Il determinante calcolato da det è una misura del fattore di scala della trasformazione lineare descritto dalla matrice. Quando il determinante è esattamente zero, la matrice è singolare e non ammette inversi.

d = det(A)

d =

     1

Alcune matrici sono quasi singolari e, nonostante il fatto che ammettano un inverso, il calcolo è soggetto a errori numerici. La funzione cond calcola il numero di condizionamento per l'inversione, che offre un'indicazione dell'accuratezza dei risultati dell'inversione della matrice. Il numero di condizionamento va da 1 per una matrice numericamente stabile a Inf per una matrice singolare.

c = cond(A)

c =

   61.9839

Raramente è necessario formare l'inverso esplicito di una matrice. Un utilizzo improprio frequente di inv si verifica quando si risolve il sistema di equazioni lineari Ax = b. Il metodo migliore per risolvere questa equazione, dal punto di vista sia del tempo di esecuzione, sia dell'accuratezza numerica, consiste nell'utilizzare l'operatore barra rovesciata x = A\b. Per maggiori informazioni consultare mldivide.

Prodotto tensoriale di Kronecker

Il prodotto di Kronecker, kron(X,Y), di due matrici è la matrice più grande formata da tutti i prodotti possibili degli elementi di X con quelli di Y. Se X è mxn e Y è pxq, allora kron(X,Y) è mpxnq. Gli elementi sono disposti in modo che ciascun elemento di X sia moltiplicato per l'intera matrice Y:

[X(1,1)*Y  X(1,2)*Y  . . .  X(1,n)*Y
                     . . .
 X(m,1)*Y  X(m,2)*Y  . . .  X(m,n)*Y]

Il prodotto di Kronecker è spesso utilizzato con matrici di zero e uno per realizzare copie ripetute di matrici più piccole. Ad esempio, se X è la matrice 2x2

X = [1   2
     3   4]

e I = eye(2,2) è la matrice di identità 2x2, allora:

kron(X,I)

ans =

     1     0     2     0
     0     1     0     2
     3     0     4     0
     0     3     0     4

kron(I,X)

ans =

     1     2     0     0
     3     4     0     0
     0     0     1     2
     0     0     3     4

A parte kron, alcune altre funzioni utili per replicare gli array sono repmat, repelem e blkdiag.

Norme di vettori e di matrici

La norma p di un vettore x,

${‖ x ‖}_{p} = {(\sum {| x_{i} |}^{p})}^{\frac{1}{p}},$

si calcola con norm(x,p). Questa operazione è definita per qualsiasi valore di p > 1, ma i valori più diffusi dip sono 1, 2 e ∞. Il valore predefinito è p = 2, che corrisponde alla lunghezza euclidea o magnitudine di un vettore:

v = [2 0 -1];
[norm(v,1) norm(v) norm(v,inf)]

ans =

    3.0000    2.2361    2.0000

La norma p di un vettore A,

${‖ A ‖}_{p} = \max_{x} \frac{{‖ A x ‖}_{p}}{{‖ x ‖}_{p}},$

può essere calcolata per p = 1, 2 e ∞ con norm(A,p). Ancora una volta il valore predefinito è p = 2:

A = pascal(3);
[norm(A,1) norm(A) norm(A,inf)]

ans =

   10.0000    7.8730   10.0000

Per calcolare la norma di ciascuna riga o colonna di una matrice utilizzare vecnorm:

vecnorm(A)

ans =

    1.7321    3.7417    6.7823

Utilizzo del calcolo multithread con le funzioni dell'algebra lineare

MATLAB supporta il calcolo multithread per diverse funzioni numeriche dell'algebra lineare e orientate agli elementi. Queste funzioni vengono eseguite automaticamente su più thread. Per una maggior rapidità di esecuzione di una funzione o espressione su più CPU, occorre soddisfare una serie di condizioni:

La funzione deve eseguire operazioni facilmente ripartibili in sezioni eseguibili contemporaneamente. Queste sezioni devono poter essere eseguite con poche comunicazioni tra i processi. Devono richiedere poche operazioni sequenziali.
Le dimensioni dei dati devono essere sufficientemente grandi, in modo che i vantaggi dell'esecuzione simultanea giustifichino il tempo richiesto per la partizione dei dati e la gestione di thread di esecuzione distinti. Ad esempio, la maggior parte delle funzioni risulta velocizzata solo quando gli array contengono almeno diverse migliaia di elementi.
L'operazione non è vincolata alla memoria; il tempo di accesso alla memoria non deve costituire la maggior parte del tempo di elaborazione. Come regola generale, le funzioni complesse sono più velocizzabili rispetto alle funzioni semplici.

Gli operatori di moltiplicazione delle matrici (X*Y) e di elevazione a potenza delle matrici (X^p) mostrano un notevole aumento della velocità sugli array con doppia precisione (nell'ordine dei 10.000 elementi). Anche le funzioni di analisi delle matrici det, rcond, hess e expm mostrano un notevole aumento della velocità sui grandi array con doppia precisione.

Argomenti complementari

Siti web esterni

Metodi matriciali di algebra lineare (Risorse didattiche MathWorks)

Sistemi di equazioni lineari

Considerazioni per i calcoli
Soluzione generale
Sistemi quadrati
Sistemi sovradeterminati
Sistemi sottodeterminati
Soluzione di sistemi con più membri sulla destra
Metodi iterativi
Calcolo multithread

Considerazioni per i calcoli

Uno dei problemi più importanti nel calcolo tecnico è la risoluzione dei sistemi di equazioni lineari simultanee.

Nella notazione matriciale, questo problema generale assume la seguente forma: Date due matrici A e b, esiste una matrice x unica, per cui Ax= b o xA= b?

È utile considerare un esempio 1x1. Ad esempio, l'equazione

7x = 21

presenta una soluzione unica?

La risposta ovviamente è sì. L'equazione presenta la soluzione unica x = 3. La soluzione viene ottenuta facilmente tramite la divisione:

x = 21/7 = 3.

La soluzione non viene ottenuta ordinariamente calcolando l'inverso di 7, cioè 7^–1= 0,142857..., e quindi moltiplicando 7^–1 per 21. Sarebbe più laborioso e, rappresentando 7^–1 con un numero finito di cifre, meno accurato. Considerazioni simili si applicano anche alle serie di equazioni lineari con più di un'incognita; MATLAB le risolve senza calcolare l'inverso della matrice.

Benché non si tratti di notazione matematica standard, MATLAB utilizza la familiare terminologia delle divisioni nel caso scalare per descrivere la soluzione di un sistema generale di equazioni simultanee. I due simboli della divisione, la barra, /, e la barra rovesciata, \, corrispondono alle due funzioni di MATLAB mrdivide e mldivide. Questi operatori sono utilizzati per le situazioni in cui la matrice ignota appare sulla sinistra o sulla destra della matrice dei coefficienti:

`x = b/A`	Denota la soluzione dell'equazione matriciale xA = b, ottenuta con `mrdivide`.
`x = A\b`	Denota la soluzione dell'equazione matriciale Ax = b, ottenuta con `mldivide`.

Si pensi alla "divisione" dei due lati dell'equazione Ax = b o xA = b per A. La matrice di coefficienti A si trova sempre nel "denominatore".

Le condizioni di compatibilità delle dimensioni per x = A\b richiedono che le due matrici A e b presentino lo stesso numero di righe. La soluzione x quindi presenta lo stesso numero di colonne di b e la dimensione delle sue righe è uguale alle dimensioni delle colonne di A. Per x = b/A i ruoli di righe e colonne risultano scambiati.

In pratica, le equazioni lineari della forma Ax = b sono più frequenti rispetto a quelle della forma xA = b. La barra rovesciata è quindi molto più utilizzata rispetto all'altro tipo di barra. La parte restante di questa sezione si concentra sull'operatore barra rovesciata; le proprietà corrispondenti dell'altro tipo di barra possono essere ricavate dall'identità:

(b/A)' = (A'\b').

La matrice dei coefficienti A non deve essere quadrata. Se A presenta le dimensioni mxn, vi sono tre casi:

m = n	Sistema quadrato. Cerca una soluzione esatta.
m > n	Sistema sovradeterminato, con più equazioni che incognite. Trova una soluzione basata sui minimi quadrati.
m < n	Sistema sottodeterminato, con meno equazioni che incognite. Trova una soluzione di base con almeno m componenti non zero.

L’algoritmo mldivide. L'operatore mldivide utilizza risolutori diversi per gestire diversi tipi di matrici dei coefficienti. Viene effettuata la diagnosi automatica dei vari casi tramite l'esame della matrice dei coefficienti. Per maggiori informazioni vedere la sezione "Algoritmi" della pagina di riferimento su mldivide.

Soluzione generale

La soluzione generale di un sistema di equazioni lineari Ax= b descrive tutte le soluzioni possibili. È possibile trovare la soluzione generale come segue:

Risolvendo il sistema omogeneo corrispondente Ax = 0. A questo scopo utilizzare il comando null, digitando null(A). Questo restituisce una base per lo spazio della soluzione verso Ax = 0. Qualsiasi soluzione è una combinazione lineare di vettori di base.
Trovando una soluzione particolare al sistema non omogeneo Ax = b.

È poi possibile scrivere qualsiasi soluzione per Ax = b come somma della soluzione particolare verso Ax =b, dal punto 2, oltre una combinazione lineare dei vettori di base dal punto 1.

La parte restante di questa sezione descrive come utilizzare MATLAB per trovare una particolare soluzione verso Ax = b, come nel punto 2.

Sistemi quadrati

La situazione più comune interessa una matrice quadrata dei coefficienti A e un singolo vettore colonna sul lato destro b.

Matrice non singolare dei coefficienti. Se la matrice A è non singolare, allora la soluzione, x = A\b, presenta le stesse dimensioni di b. Ad esempio:

A = pascal(3);
u = [3; 1; 4];
x = A\u

x =
      10
     -12
       5

È possibile confermare che A*x è esattamente pari a u.

Se A e b sono quadrate e delle stesse dimensioni, anche x= A\b presenta le stesse dimensioni:

b = magic(3);
X = A\b

X =
      19    -3    -1
     -17     4    13
       6     0    -6

È possibile confermare che A*x è esattamente uguale a b.

Entrambi gli esempi generano soluzioni esatte e composte da numeri interi. Questo accade perché si è scelta una matrice dei coefficienti pascal(3), che è una matrice a rank pieno (non singolare).

Matrice singolare dei coefficienti. Una matrice quadrata A è singolare se non presenta colonne linearmente indipendenti. Se A è singolare, la soluzione di Ax = b non esiste oppure non è unica. L'operatore barra rovesciata, A\b, genera un avviso se A è quasi singolare o se rileva una singolarità esatta.

Se A è singolare e Ax = b ha una soluzione, è possibile trovare una soluzione particolare non unica digitando

P = pinv(A)*b

pinv(A) è una pseudo-inversa di A. Se Ax = b non presenta una soluzione esatta, pinv(A) restituisce una soluzione basata sui minimi quadrati.

Ad esempio:

A = [ 1     3     7
     -1     4     4
      1    10    18 ]

è singolare, come è possibile verificare digitando

rank(A)

ans =

     2

Poiché A non è a rank pieno, presenta alcuni singoli valori uguali a zero.

Soluzioni esatte. Per b =[5;2;12], l'equazione Ax = b presenta una soluzione esatta, data da

pinv(A)*b

ans =
    0.3850
   -0.1103
    0.7066

Verificare che pinv(A)*b sia una soluzione esatta digitando

A*pinv(A)*b

ans =
    5.0000
    2.0000
   12.0000

Soluzioni basate sui minimi quadrati. Se tuttavia b = [3;6;0], Ax = b non presenta una soluzione esatta. In questo caso pinv(A)*b restituisce una soluzione basata sui minimi quadrati. Digitando

A*pinv(A)*b

ans =
   -1.0000
    4.0000
    2.0000

non si ottiene nuovamente il vettore originale b.

È possibile determinare se Ax = b presenta una soluzione esatta trovando la forma a gradini ridotta della riga della matrice aumentata [A b]. Per eseguire questa operazione sull'esempio, digitare

rref([A b])
ans =
    1.0000         0    2.2857         0
         0    1.0000    1.5714         0
         0         0         0    1.0000

Poiché la riga inferiore contiene tutti zeri, a eccezione dell'ultimo elemento, l'equazione non ha soluzione. In questo caso pinv(A) restituisce una soluzione basata sui minimi quadrati.

Sistemi sovradeterminati

Apri live script

Questo esempio dimostra quanto siano frequenti i sistemi sovradeterminati in diversi tipi di adattamento delle curve a dati sperimentali.

Una quantità y viene misurata in diversi valori di tempo t per produrre le osservazioni sotto riportate. Le seguenti dichiarazioni consentono di inserire i dati e di visualizzarli in una tabella.

t = [0 .3 .8 1.1 1.6 2.3]';
y = [.82 .72 .63 .60 .55 .50]';
B = table(t,y)

B=6×2 table
     t      y  
    ___    ____

      0    0.82
    0.3    0.72
    0.8    0.63
    1.1     0.6
    1.6    0.55
    2.3     0.5

Provare a modellare i dati con una funzione di decadimento esponenziale

$y (t) = c_{1} + c_{2} e^{- t}$ .

L'equazione precedente dichiara che il vettore y dovrebbe essere approssimato da una combinazione lineare di altri due vettori. Uno è un vettore costante contenente tutti uno, l'altro è il vettore con i componenti exp(-t). È possibile calcolare i coefficienti incogniti $c_{1}$ e $c_{2}$ con un adattamento basato sui minimi quadrati, che minimizza la somma dei quadrati delle deviazioni dei dati dal modello. Vi sono sei equazioni in due incognite, rappresentate da una matrice 6x2.

E = [ones(size(t)) exp(-t)]

E = 6×2

    1.0000    1.0000
    1.0000    0.7408
    1.0000    0.4493
    1.0000    0.3329
    1.0000    0.2019
    1.0000    0.1003

Utilizzare l'operatore barra rovesciata per ottenere la soluzione basata sui minimi quadrati.

c = E\y

In altri termini, l'adattamento basato sui minimi quadrati dei dati è

$y (t) = 0.4760 + 0.3413 e^{- t} .$

Le dichiarazioni seguenti valutano il modello in corrispondenza di incrementi con spaziamento regolare in t, quindi tracciano il risultato insieme ai dati originali:

T = (0:0.1:2.5)';
Y = [ones(size(T)) exp(-T)]*c;
plot(T,Y,'-',t,y,'o')

Figure contains an axes object. The axes object contains 2 objects of type line. One or more of the lines displays its values using only markers

E*c non è esattamente uguale a y, ma la differenza potrebbe risultare inferiore rispetto agli errori di misurazione nei dati originali.

Una matrice rettangolare A presenta una carenza di rank se non ha colonne linearmente indipendenti. Se A presenta una carenza di rank, la soluzione basata sui minimi quadrati di AX = B non è unica. A\B genera un errore se A presenta una carenza di rank e genera una soluzione basata sui minimi quadrati. È possibile utilizzare lsqminnorm per trovare la soluzione di X con norma minima tra tutte le soluzioni.

Sistemi sottodeterminati

Questo esempio mostra come la soluzione dei sistemi sottodeterminati non sia unica. I sistemi lineari sottodeterminati presentano più incognite che equazioni. L'operazione di divisione sulla sinistra della matrice in MATLAB trova una soluzione di base basata sui minimi quadrati, che ha al massimo m componenti non zero per una matrice dei coefficienti mxn.

Ecco un esempio casuale:

R = [6 8 7 3; 3 5 4 1]
rng(0);
b = randi(8,2,1)

R =

       6              8              7              3       
       3              5              4              1       


b =

       7       
       8

Il sistema lineare Rp = b implica due equazioni in quattro incognite. Poiché la matrice dei coefficienti contiene numeri interi di valore ridotto, è appropriato usare il comando format per visualizzare la soluzione in formato razionale. Questa particolare soluzione si ottiene con

format rat
p = R\b

Uno dei componenti non zero è p(2), perché R(:,2) è la colonna di R con la norma più grande. L'altro componente non zero è p(4), perché R(:,4) risulta dominante dopo l'eliminazione di R(:,2).

La soluzione generale completa del sistema sottodeterminato può essere caratterizzata aggiungendo p a una combinazione lineare arbitraria dei vettori spazio nulli, individuabili tramite la funzione null con un'opzione che richiede una base razionale.

Z = null(R,'r')

Z =

      -1/2           -7/6     
      -1/2            1/2     
       1              0       
       0              1

È possibile confermare che R*Z è zero e che R*x - b residuo è piccolo per qualsiasi vettore x, dove

x = p + Z*q

Poiché le colonne di Z sono i vettori spazio nulli, il prodotto Z*q è una combinazione lineare di questi vettori:

$Z q = (\begin{matrix} {\overset{⇀}{x}}_{1} & {\overset{⇀}{x}}_{2} \end{matrix}) (\begin{matrix} u \\ w \end{matrix}) = u {\overset{⇀}{x}}_{1} + w {\overset{⇀}{x}}_{2} .$

Come esempio illustrativo, scegliere un q arbitrario e costruire x.

q = [-2; 1];
x = p + Z*q;

Calcolare la norma del residuo.

format short
norm(R*x - b)

ans =

   2.6645e-15

Quando vi sono molte soluzioni disponibili, quella con la norma minima risulta di particolare interesse. È possibile usare lsqminnorm per calcolare la soluzione con norma minima basata sui minimi quadrati. Questa soluzione presenta il valore minimo possibile per norm(p).

p = lsqminnorm(R,b)

Soluzione di sistemi con più membri sulla destra

Alcuni problemi riguardano la soluzione di sistemi lineari che presentano la stessa matrice dei coefficienti A, ma diversi membri sulla destra b. Quando sono disponibili contemporaneamente diversi valori di b, è possibile costruire b come matrice con diverse colonne e risolvere contemporaneamente tutti i sistemi di equazioni utilizzando un solo comando barra rovesciata: X = A\[b1 b2 b3 …].

A volte, tuttavia, i diversi valori di b non sono tutti disponibili contemporaneamente, il che significa che è necessario risolvere di seguito diversi sistemi di equazioni. Quando si risolve uno di questi sistemi di equazioni con la barra (/) o la barra rovesciata (\), l'operatore fattorizza la matrice dei coefficienti A e utilizza questa decomposizione in matrice per calcolare la soluzione. A ogni successiva risoluzione di un sistema di equazioni simili con diverso b, tuttavia, l'operatore calcola la stessa decomposizione di A, che è un calcolo ridondante.

La soluzione a questo problema consiste nel calcolare preventivamente la decomposizione di A, quindi riutilizzare i fattori per risolvere i diversi valori di b. In pratica, però, il calcolo preventivo della decomposizione in questo modo può risultare difficile, perché è necessario sapere quale decomposizione calcolare (LU, LDL, Cholesky e così via) e anche come moltiplicare i fattori per risolvere il problema. Ad esempio, con la decomposizione LU è necessario risolvere due sistemi lineari per risolvere il sistema originale Ax = b:

[L,U] = lu(A);
x = U \ (L \ b);

Il metodo consigliato per la risoluzione di sistemi lineari con più membri sulla destra consecutivi, invece, consiste nell'utilizzare gli oggetti decomposition. Questi oggetti consentono di sfruttare i vantaggi delle performance del calcolo preventivo della decomposizione in matrice, ma non richiedono di sapere come utilizzare i fattori delle matrici. È possibile sostituire la precedente decomposizione LU con:

dA = decomposition(A,'lu');
x = dA\b;

Se non si è certi del tipo di decomposizione da utilizzare, decomposition(A) sceglie il tipo corretto in base alle proprietà di A, in modo simile a quanto viene effettuato con l'operatore barra retroversa.

Ecco un semplice test dei possibili vantaggi di questo approccio in termini di performance. Il test risolve 100 volte lo stesso sistema lineare sparso utilizzando sia (\) che decomposition.

n = 1e3;
A = sprand(n,n,0.2) + speye(n);
b = ones(n,1);

% Backslash solution
tic
for k = 1:100
    x = A\b;
end
toc

Elapsed time is 9.006156 seconds.

% decomposition solution
tic
dA = decomposition(A);
for k = 1:100
    x = dA\b;
end
toc

Elapsed time is 0.374347 seconds.

Per questo problema, la soluzione decomposition è molto più rapida rispetto all'uso della sola barra retroversa, pur con una sintassi semplice.

Metodi iterativi

Se la matrice dei coefficienti A è grande e sparsa, solitamente i metodi di fattorizzazione non sono efficienti. I metodi iterativi generano una serie di soluzioni approssimative. MATLAB offre diversi metodi iterativi per la gestione delle matrici di input grandi e sparse.

Funzione	Descrizione
`pcg`	Metodo del gradiente coniugato precondizionato. Questo metodo è appropriato per una matrice A con coefficiente hermitiano definito positivo.
`bicg`	Metodo del gradiente biconiugato
`bicgstab`	Metodo stabilizzato del gradiente biconiugato
`bicgstabl`	Metodo BiCGStab(l)
`cgs`	Metodo quadrato gradiente coniugato
`gmres`	Metodo del residuo minimo generalizzato
`lsqr`	Metodo LSQR
`minres`	Metodo del residuo minimo. Questo metodo è appropriato per una matrice A con coefficiente hermitiano.
`qmr`	Metodo del residuo quasi-minimo
`symmlq`	Metodo LQ simmetrico
`tfqmr`	Metodo QMR senza trasposizione

Calcolo multithread

La funzione deve eseguire operazioni facilmente ripartibili in sezioni eseguibili contemporaneamente. Queste sezioni devono poter essere eseguite con poche comunicazioni tra i processi. Devono richiedere poche operazioni sequenziali.
Le dimensioni dei dati devono essere sufficientemente grandi, in modo che i vantaggi dell'esecuzione simultanea giustifichino il tempo richiesto per la partizione dei dati e la gestione di thread di esecuzione distinti. Ad esempio, la maggior parte delle funzioni risulta velocizzata solo quando gli array contengono almeno diverse migliaia di elementi.
L'operazione non è vincolata alla memoria; il tempo di accesso alla memoria non deve costituire la maggior parte del tempo di elaborazione. Come regola generale, le funzioni complesse sono più velocizzabili rispetto alle funzioni semplici.

inv, lscov, linsolve e mldivide mostrano un notevole aumento della velocità sugli array a doppia precisione (nell'ordine di almeno 10.000 elementi) quando la modalità multithread è attiva.

Vedi anche

mldivide | mrdivide | pinv | decomposition | lsqminnorm

Argomenti complementari

Siti web esterni

Metodi matriciali di algebra lineare (Risorse didattiche MathWorks)

Fattorizzazioni

Introduzione
Fattorizzazione di Cholesky
Fattorizzazione LU
Fattorizzazione QR
Uso del calcolo multithread per la fattorizzazione

Introduzione

Tutte e tre le fattorizzazioni di matrici discusse in questa sezione utilizzano le matrici triangolari, dove tutti gli elementi sopra o sotto la diagonale sono pari a zero. I sistemi di equazioni lineari che interessano matrici triangolari possono essere risolti facilmente e semplicemente utilizzando la sostituzione in avanti o all'indietro.

Fattorizzazione di Cholesky

La fattorizzazione di Cholesky esprime una matrice simmetrica come prodotto di una matrice triangolare e il suo trasposto

A = R′R,

dove R è una matrice triangolare superiore.

Non tutte le matrici simmetriche possono essere fattorizzate in questo modo; le matrici che presentano questa fattorizzazione sono denominate come positive definite. Questo implica che tutti gli elementi della diagonale di A siano positivi e che gli elementi della contro-diagonale "non siano troppo grandi". Le matrici di Pascal rappresentano un esempio interessante. In tutto questo capitolo si è utilizzata una matrice Pascal 3x3 di esempio A. Passare momentaneamente a 6x6:

A = pascal(6)

A =
       1     1     1     1     1     1
       1     2     3     4     5     6
       1     3     6    10    15    21
       1     4    10    20    35    56
       1     5    15    35    70   126
       1     6    21    56   126   252

Gli elementi di A sono coefficienti binomiali. Ciascun elemento è la somma degli elementi vicini a nord e ovest. La fattorizzazione di Cholesky è

R = chol(A)

R =
     1     1     1     1     1     1
     0     1     2     3     4     5
     0     0     1     3     6    10
     0     0     0     1     4    10
     0     0     0     0     1     5
     0     0     0     0     0     1

Gli elementi sono di nuovo coefficienti binomiali. Il fatto che R'*R sia uguale a A dimostra un'identità che interessa somme di prodotti di coefficienti binomiali.

Nota

La fattorizzazione di Cholesky si applica anche alle matrici complesse. Qualsiasi matrice complessa con fattorizzazione di Cholesky soddisfa

A′ = A

ed è denominata come definita positiva hermitiana.

La fattorizzazione di Cholesky consente di sostituire il sistema lineare

Ax = b

con

R′Rx = b.

Dato che l'operatore barra rovesciata riconosce i sistemi triangolari, questo può essere risolto rapidamente nell'ambiente di MATLAB con

x = R\(R'\b)

Se A è nxn, la complessità di calcolo di chol(A) è O(n³), ma la complessità delle successive soluzioni con barra retroversa è solo O(n²).

Fattorizzazione LU

La fattorizzazione LU, o eliminazione gaussiana, esprime una qualsiasi matrice quadrata A come prodotto di una permutazione di una matrice triangolare inferiore e una matrice triangolare superiore

A = LU,

dove L è una permutazione di una matrice triangolare inferiore con elementi uno sulla sua diagonale e U è una matrice triangolare superiore.

Le permutazioni sono necessarie per motivi sia teorici che di calcolo. Non è possibile esprimere la matrice

$[\begin{matrix} 0 & 1 \\ 1 & 0 \end{matrix}]$

come prodotto di matrici triangolari senza scambiarne le due righe. Benché la matrice

$[\begin{matrix} ε & 1 \\ 1 & 0 \end{matrix}]$

sia esprimibile come prodotto di matrici triangolari, quando ε è piccolo, gli elementi nei fattori sono grandi e amplificano gli errori; quindi, benché non strettamente necessarie, le permutazioni sono consigliabili. La pivotazione parziale garantisce che gli elementi di L siano limitati da uno in magnitudine e che gli elementi di U non siano molto più grandi di quelli in A.

Ad esempio:

[L,U] = lu(B)

L =
    1.0000         0         0
    0.3750    0.5441    1.0000
    0.5000    1.0000         0

U =
    8.0000    1.0000    6.0000
         0    8.5000   -1.0000
         0         0    5.2941

La fattorizzazione LU di A consente di risolvere rapidamente il sistema lineare

A*x = b

con

x = U\(L\b)

Determinanti e inverse sono calcolate dalla fattorizzazione LU con

det(A) = det(L)*det(U)

inv(A) = inv(U)*inv(L)

È anche possibile calcolare le determinanti con det(A) = prod(diag(U)), anche se i segni delle determinanti possono essere invertiti.

Fattorizzazione QR

Una matrice ortogonale, o una matrice con colonne ortonormali, è una matrice reale con colonne tutte di lunghezza unitaria e perpendicolari l'una con l'altra. Se Q è ortogonale, allora

Q^TQ = I,

dove I è la matrice di identità.

Le matrici ortogonali più semplici sono rotazioni di coordinate bidimensionali:

$[\begin{matrix} \cos (θ) & \sin (θ) \\ - \sin (θ) & \cos (θ) \end{matrix}] .$

Per le matrici complesse, il termine corrispondente è unitaria. Le matrici ortogonali e unitarie sono utili nel calcolo numerico perché conservano la lunghezza, conservano gli angoli e non amplificano gli errori.

La fattorizzazione ortogonale, o QR, esprime una qualsiasi matrice rettangolare come prodotto di una matrice ortogonale o unitaria e di una matrice triangolare superiore. È inoltre possibile effettuare una permutazione di colonna:

A = QR

oppure

AP = QR,

dove Q è ortogonale o unitario, R è un triangolo superiore e P è una permutazione.

Vi sono quattro varianti della fattorizzazione QR: complete o economy size e con o senza permutazioni di colonne.

I sistemi lineari sovradeterminati interessano una matrice rettangolare con più righe che colonne, cioè mxn con m > n. La fattorizzazione QR a piene dimensioni produce un quadrato ortogonale Q, mxm e un triangolo superiore rettangolare R,mxn:

C=gallery('uniformdata',[5 4], 0);
[Q,R] = qr(C)

Q =

    0.6191    0.1406   -0.1899   -0.5058    0.5522
    0.1506    0.4084    0.5034    0.5974    0.4475
    0.3954   -0.5564    0.6869   -0.1478   -0.2008
    0.3167    0.6676    0.1351   -0.1729   -0.6370
    0.5808   -0.2410   -0.4695    0.5792   -0.2207


R =

    1.5346    1.0663    1.2010    1.4036
         0    0.7245    0.3474   -0.0126
         0         0    0.9320    0.6596
         0         0         0    0.6648
         0         0         0         0

In molti casi, le ultime colonne m – n di Q non sono necessarie perché vengono moltiplicate per gli zeri nella parte inferiore di R. La fattorizzazione QR economy size produce un rettangolo Q mxn con colonne ortonormali e un quadrato triangolare superiore R nxn. Per l'esempio 5x4 non si tratta di un vantaggio importante, mentre per le matrici più grandi e molto rettangolari il risparmio in termini di tempo e memoria può essere decisamente notevole:

[Q,R] = qr(C,0)	
Q =

    0.6191    0.1406   -0.1899   -0.5058
    0.1506    0.4084    0.5034    0.5974
    0.3954   -0.5564    0.6869   -0.1478
    0.3167    0.6676    0.1351   -0.1729
    0.5808   -0.2410   -0.4695    0.5792


R =

    1.5346    1.0663    1.2010    1.4036
         0    0.7245    0.3474   -0.0126
         0         0    0.9320    0.6596
         0         0         0    0.6648

Contrariamente alla fattorizzazione LU, la fattorizzazione QR non richiede pivotazione o permutazioni. Ma una permutazione opzionale delle colonne, attivata dalla presenza di un terzo argomento di output, è utile per la rilevazione di singolarità o di carenze di rango. In qualsiasi passaggio della fattorizzazione, viene usata come base la colonna delle restanti matrici non sottoposte a fattorizzazione con la norma più grande. Questo garantisce che gli elementi diagonali di R si presentino in ordine decrescente e che eventuali dipendenze lineari tra le colonne vengano quasi certamente rivelate da un esame di questi elementi. Per il breve esempio qui riportato, la seconda colonna di C presenta una norma maggiore della prima, quindi le due colonne vengono scambiate:

[Q,R,P] = qr(C)

Q =
   -0.3522    0.8398   -0.4131
   -0.7044   -0.5285   -0.4739
   -0.6163    0.1241    0.7777

R =
  -11.3578   -8.2762
         0    7.2460
         0         0

P =
     0     1
     1     0

Quando le permutazioni economy size e di colonna vengono combinate, il terzo argomento di output è un vettore di permutazione, non una matrice di permutazione:

[Q,R,p] = qr(C,0)

Q =
   -0.3522    0.8398
   -0.7044   -0.5285
   -0.6163    0.1241

R =
  -11.3578   -8.2762
         0    7.2460


p =
     2     1

La fattorizzazione QR trasforma un sistema lineare sovradeterminato in un sistema triangolare equivalente. L'espressione

norm(A*x - b)

è uguale a

norm(Q*R*x - b)

La moltiplicazione con matrici ortogonali osserva la norma euclidea, quindi questa espressione è anche uguale a

norm(R*x - y)

dove y = Q'*b. Poiché le ultime righe m-n di R sono degli zeri, questa espressione si suddivide in due parti:

norm(R(1:n,1:n)*x - y(1:n))

norm(y(n+1:m))

Quando A presenta un rank pieno, è possibile risolverlo per x in modo che la prima di queste espressioni sia zero. La seconda espressione dà poi la norma del residuo. Quando A non presenta un rango completo, la struttura triangolare di R consente di trovare una soluzione di base per il problema con minimi quadrati.

Uso del calcolo multithread per la fattorizzazione

Il software MATLAB supporta il calcolo multithread per diverse funzioni numeriche dell'algebra lineare e orientate agli elementi. Queste funzioni vengono eseguite automaticamente su più thread. Per una maggior rapidità di esecuzione di una funzione o espressione su più CPU, occorre soddisfare una serie di condizioni:

La funzione deve eseguire operazioni facilmente ripartibili in sezioni eseguibili contemporaneamente. Queste sezioni devono poter essere eseguite con poche comunicazioni tra i processi. Devono richiedere poche operazioni sequenziali.
Le dimensioni dei dati devono essere sufficientemente grandi, in modo che i vantaggi dell'esecuzione simultanea giustifichino il tempo richiesto per la partizione dei dati e la gestione di thread di esecuzione distinti. Ad esempio, la maggior parte delle funzioni risulta velocizzata solo quando gli array contengono almeno diverse migliaia di elementi.
L'operazione non è vincolata alla memoria; il tempo di accesso alla memoria non deve costituire la maggior parte del tempo di elaborazione. Come regola generale, le funzioni complesse sono più velocizzabili rispetto alle funzioni semplici.

lu e qr mostrano un notevole aumento della velocità sugli array a doppia precisione (nell'ordine di 10.000 elementi).

Vedi anche

lu | qr | chol

Argomenti complementari

Le matrici nell'ambiente di MATLAB

Potenze ed esponenti

Apri live script

Questo argomento mostra come calcolare potenze ed esponenti di matrici con diversi metodi.

Potenze di numeri interi positivi

Se A è una matrice quadrata e p è un numero positivo intero, allora A^p moltiplica effettivamente A per se stesso p-1 volte. Ad esempio:

ans = 3×3

     3     6    10
     6    14    25
    10    25    46

Potenze inverse e frazionali

Se A è una matrice quadrata e non singolare, allora A^(-p) moltiplica effettivamente inv(A) per se stesso p-1 volte.

A^(-3)

ans = 3×3

  145.0000 -207.0000   81.0000
 -207.0000  298.0000 -117.0000
   81.0000 -117.0000   46.0000

MATLAB® calcola inv(A) e A^(-1) con lo stesso algoritmo, quindi i risultati sono esattamente identici. Sia inv(A) che A^(-1) generano avvisi se la matrice è vicina a essere singolare.

isequal(inv(A),A^(-1))

ans = logical
   1

È consentito anche l'utilizzo di potenze frazionarie, come A^(2/3). I risultati ottenuti con potenze frazionarie dipendono dalla distribuzione degli autovalori della matrice.

A^(2/3)

ans = 3×3

    0.8901    0.5882    0.3684
    0.5882    1.2035    1.3799
    0.3684    1.3799    3.1167

Potenze elemento per elemento

L'operatore .^ calcola le potenze elemento per elemento. Ad esempio, per calcolare il quadrato di ogni elemento in una matrice è possibile utilizzare A.^2.

A.^2

ans = 3×3

     1     1     1
     1     4     9
     1     9    36

Radici quadrate

La funzione sqrt rappresenta un metodo utile per calcolare la radice quadrata di ogni elemento in una matrice. Un altro modo per eseguire questo calcolo consiste nell'utilizzare A.^(1/2).

sqrt(A)

ans = 3×3

    1.0000    1.0000    1.0000
    1.0000    1.4142    1.7321
    1.0000    1.7321    2.4495

Per altre radici è possibile utilizzare nthroot. Ad esempio, calcolare A.^(1/3).

nthroot(A,3)

ans = 3×3

    1.0000    1.0000    1.0000
    1.0000    1.2599    1.4422
    1.0000    1.4422    1.8171

Queste radici orientate agli elementi sono diverse dalle radici quadrate di una matrice, che calcolano una seconda matrice $B$ come $A = BB$ . La funzione sqrtm(A) calcola A^(1/2) con un algoritmo più accurato. La m in sqrtm distingue questa funzione da sqrt(A) che, come A.^(1/2), utilizza un metodo 'elemento per elemento'.

B = sqrtm(A)

B = 3×3

    0.8775    0.4387    0.1937
    0.4387    1.0099    0.8874
    0.1937    0.8874    2.2749

B^2

ans = 3×3

    1.0000    1.0000    1.0000
    1.0000    2.0000    3.0000
    1.0000    3.0000    6.0000

Basi scalari

Oltre a elevare una matrice a potenza, è anche possibile elevare uno scalare alla potenza di una matrice.

2^A

ans = 3×3

   10.4630   21.6602   38.5862
   21.6602   53.2807   94.6010
   38.5862   94.6010  173.7734

Quando si eleva uno scalare alla potenza di una matrice, MATLAB utilizza gli autovalori e gli autovettori della matrice per calcolare la potenza della matrice stessa. Se [V,D] = eig(A), allora $2^{A} = {V 2}^{D} V^{- 1}$ .

[V,D] = eig(A);
V*2^D*V^(-1)

ans = 3×3

   10.4630   21.6602   38.5862
   21.6602   53.2807   94.6010
   38.5862   94.6010  173.7734

Esponenziali di matrice

Gli esponenziali di matrici rappresentano un caso speciale per elevare uno scalare a una potenza di matrice. La base di un esponenziale di matrice è il numero di Eulero e = exp(1).

e = exp(1);
e^A

ans = 3×3
10³ ×

    0.1008    0.2407    0.4368
    0.2407    0.5867    1.0654
    0.4368    1.0654    1.9418

La funzione expm rappresenta un metodo più utile per calcolare gli esponenziali di matrice.

expm(A)

ans = 3×3
10³ ×

    0.1008    0.2407    0.4368
    0.2407    0.5867    1.0654
    0.4368    1.0654    1.9418

Gli esponenziali di matrice possono essere calcolati in diversi modi. Per maggiori informazioni vedere Matrix Exponentials.

Operazioni con numeri piccoli

Le funzioni log1p e expm1 di MATLAB calcolano accuratamente $\log (1 + x)$ e $e^{x} - 1$ per valori molto piccoli di $x$ . Ad esempio, se si prova ad aggiungere un numero inferiore rispetto alla precisione di macchina di 1, il risultato viene arrotondato a 1.

log(1+eps/2)

ans = 0

log1p è tuttavia in grado di restituire una risposta più accurata.

log1p(eps/2)

ans = 1.1102e-16

Come accade per $e^{x} - 1$ , se è molto piccolo, $x$ viene arrotondato a zero.

exp(eps/2)-1

ans = 0

Ancora una volta, expm1 è in grado di restituire una risposta più accurata.

expm1(eps/2)

ans = 1.1102e-16

Vedi anche

Argomenti complementari

Decomposizione di autovalori

Un autovalore e un autovettore di una matrice quadrata A sono, rispettivamente, uno scalare λ e un vettore υ non zero che soddisfano

Aυ = λυ.

Con gli autovalori sulla diagonale di una matrice di diagonale Λ e i corrispettivi autovettori che formano le colonne di una matrice V, si ottiene

AV = VΛ.

Se V è non singolare, questo diventa la decomposizione di autovalori

A = VΛV^–1.

Un buon esempio è rappresentato dalla matrice di coefficienti dell'equazione differenziale dx/dt = Ax:

A =
     0    -6    -1
     6     2   -16
    -5    20   -10

La soluzione a questa equazione viene espressa in termini della matrice esponenziale x(t) = e^tAx(0). La dichiarazione

lambda = eig(A)

genera un vettore di colonna contenente l'autovalore A. Per questa matrice, gli autovalori sono complessi:

lambda =
     -3.0710         
     -2.4645+17.6008i
     -2.4645-17.6008i

La parte reale di ciascuno degli autovalori è negativa, quindi e^λt si avvicina a zero all'aumentare di t. La parte immaginaria non zero di due degli autovalori, ±ω, contribuisce alla componente oscillatoria, sin(ωt), della soluzione dell'equazione differenziale.

Con due argomenti di output, eig calcola gli autovettori e memorizza gli autovalori in una matrice di diagonali:

[V,D] = eig(A)

V =
  -0.8326         0.2003 - 0.1394i   0.2003 + 0.1394i
  -0.3553        -0.2110 - 0.6447i  -0.2110 + 0.6447i
  -0.4248        -0.6930            -0.6930          

D =
  -3.0710                 0                 0         
        0           -2.4645+17.6008i        0         
        0                 0           -2.4645-17.6008i

Il primo autovettore è reale, mentre gli altri due vettori sono coniugati complessi reciproci. Tutti e tre i vettori sono normalizzati per presentare una lunghezza euclidea, norm(v,2), pari a uno.

La matrice V*D*inv(V), che può essere scritta più brevemente come V*D/V, rientra nell'errore di arrotondamento di A. inv(V)*A*V, o V\A*V, rientra nell'errore di arrotondamento di D.

Autovalori multipli

Alcune matrici non dispongono della decomposizione degli autovettori. Si tratta di matrici non diagonizzabili. Ad esempio:

A = [ 1    -2    1 
      0     1    4 
      0     0    3 ]

Per questa matrice

[V,D] = eig(A)

genera

V =

    1.0000    1.0000   -0.5571
         0    0.0000    0.7428
         0         0    0.3714


D =

     1     0     0
     0     1     0
     0     0     3

Vi è un doppio autovalore in λ = 1. La prima e la seconda colonna di V coincidono. Per questa matrice non esiste una serie completa di autovettori linearmente indipendenti.

Decomposizione di Schur

Molti calcoli avanzati sulle matrici non richiedono la decomposizione degli autovalori. Si basano invece sulla decomposizione di Schur

A = USU ′ ,

in cui U è una matrice ortogonale e S è una matrice triangolare superiore di blocchi con blocchi 1x1 e 2x2 sulla diagonale. Gli autovalori sono rivelati dagli elementi della diagonale e dai blocchi di S, mentre le colonne di U forniscono una base ortogonale, che presenta proprietà numeriche assai superiori rispetto a una serie di autovettori.

Ad esempio, confrontare l'autovalore e le decomposizioni di Schur di questa matrice difettosa:

A = [ 6    12    19 
     -9   -20   -33 
      4     9    15 ];

[V,D] = eig(A)

V =

  -0.4741 + 0.0000i  -0.4082 - 0.0000i  -0.4082 + 0.0000i
   0.8127 + 0.0000i   0.8165 + 0.0000i   0.8165 + 0.0000i
  -0.3386 + 0.0000i  -0.4082 + 0.0000i  -0.4082 - 0.0000i


D =

  -1.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   1.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   1.0000 - 0.0000i

[U,S] = schur(A)

U =

   -0.4741    0.6648    0.5774
    0.8127    0.0782    0.5774
   -0.3386   -0.7430    0.5774


S =

   -1.0000   20.7846  -44.6948
         0    1.0000   -0.6096
         0    0.0000    1.0000

La matrice A è difettosa perché non dispone di una serie completa di autovettori linearmente indipendenti (la seconda e la terza colonna di V coincidono). Poiché non tutte le colonne di V sono linearmente indipendenti, vi è un grande numero di condizioni di circa ~1e8. schur è tuttavia in grado di calcolare tre diversi vettori di base in U. Poiché U è ortogonale, allora cond(U) = 1.

La matrice S presenta l'autovalore reale come primo elemento della diagonale e l'autovalore ripetuto rappresentato dal blocco 2x2 inferiore a destra. Anche gli autovalori del blocco 2x2 sono autovalori di A:

eig(S(2:3,2:3))

ans =

   1.0000 + 0.0000i
   1.0000 - 0.0000i

Vedi anche

eig | schur

Argomenti complementari

Siti web esterni

Metodi matriciali di algebra lineare (Risorse didattiche MathWorks)

Valori singolari

Un valore singolare e i vettori singolari corrispondenti di una matrice rettangolare A sono rispettivamente un σ scalare e una coppia di vettori u e v che soddisfano

$\begin{array}{l} A v = σ u \\ A^{H} u = σ v, \end{array}$

in cui $A^{H}$ è il trasposto hermitiano di A. I vettori singolari u e v solitamente sono scalati per presentare una norma di 1. Inoltre, se u e v sono vettori singolari di A, allora anche -u e -v sono vettori singolari di A.

I valori singolari di σ sono sempre reali e non negativi, anche se A è complesso. Con i valori singolari in una matrice di diagonali Σ e i corrispondenti vettori singolari che formano le colonne di due matrici ortogonali U e V, si ottengono le equazioni

$\begin{array}{l} A V = U Σ \\ A^{H} U = V Σ . \end{array}$

Poiché U e V sono matrici unitarie, moltiplicando la prima equazione per $V^{H}$ sulla destra si ottiene l'equazione di decomposizione in valori singolari

$A = U Σ V^{H} .$

La decomposizione completa in valori singolari di una matrice mxn comporta:

Matrice U mxm
Matrice Σ mxn
Matrice Vnxn

In altre parole, U e V sono entrambi quadrate e Σ presenta le stesse dimensioni di A. Se A presenta molte più righe che colonne (m > n), allora la matrice U mxm che ne risulta è grande. La maggior parte delle colonne in U, tuttavia, è moltiplicata per zero in Σ. In questo caso, la decomposizione economy sized consente di risparmiare tempo e spazio di memorizzazione producendo una U mxn, una Σ nxn e la stessa V:

In the economy-sized decomposition, columns in U can be ignored if they multiply zeros in the diagonal matrix of singular values.

La decomposizione in autovalori è lo strumento appropriato per analizzare una matrice quando quest'ultima rappresenta una mappatura da uno spazio di vettori in se stessa, come avviene per una normale equazione differenziale. La decomposizione in valori singoli è invece lo strumento appropriato per analizzare una mappatura da uno spazio di vettori in un altro spazio di vettori, possibilmente di dimensioni diverse. Molti sistemi di equazioni lineari simultanee rientrano in questa seconda categoria.

Se A è singolare, simmetrica e definita positiva, allora le sue decomposizioni in autovalori e in valori singolari coincidono. Poiché A, invece, si discosta dalla simmetria e dalla definitività positiva, la differenza tra le due decomposizioni aumenta. In particolare, la decomposizione in valori singolari di una matrice reale è sempre reale, mentre la decomposizione in autovalori di una matrice reale non simmetrica potrebbe essere complessa.

Per la matrice dell'esempio

A = [9     4
     6     8
     2     7];

la decomposizione completa in valori singolari è

[U,S,V] = svd(A)

U =

   -0.6105    0.7174    0.3355
   -0.6646   -0.2336   -0.7098
   -0.4308   -0.6563    0.6194


S =

   14.9359         0
         0    5.1883
         0         0


V =

   -0.6925    0.7214
   -0.7214   -0.6925

È possibile verificare che U*S*V' è pari a A entro i margini di arrotondamento degli errori. Per questo problema di piccola entità, la decomposizione economy size è solo leggermente più piccola.

[U,S,V] = svd(A,"econ")

U =

   -0.6105    0.7174
   -0.6646   -0.2336
   -0.4308   -0.6563


S =

   14.9359         0
         0    5.1883


V =

   -0.6925    0.7214
   -0.7214   -0.6925

Ancora una volta, U*S*V' è uguale a A entro i margini di arrotondamento degli errori.

Calcolo della SVD in batch

Qualora si debba decomporre un grande gruppo di matrici della stessa dimensione, eseguire tutte le decomposizioni in un loop con svd risulta essere inefficiente. È invece possibile concatenare tutte le matrici in un array multidimensionale e utilizzare pagesvd per eseguire le decomposizioni ai valori singolari su tutte le pagine dell’array, tramite un'unica chiamata di funzione.

Funzione	Utilizzo
`pagesvd`	Utilizzare `pagesvd` per eseguire le decomposizioni ai valori singolari sulle pagine di un array multidimensionale. Si tratta di un modo efficiente per eseguire la SVD su un grande gruppo di matrici della stessa dimensione.

Ad esempio, si consideri un gruppo di tre matrici 2 per 2. Concatenare le matrici in un array di 2 per 2 per 3 utilizzando la funzione cat.

A = [0 -1; 1 0];
B = [-1 0; 0 -1];
C = [0 1; -1 0];
X = cat(3,A,B,C);

Ora, utilizzare pagesvd per eseguire le tre decomposizioni contemporaneamente.

[U,S,V] = pagesvd(X);

Per ciascuna pagina di X, sono presenti le pagine corrispondenti negli output U, S e V. Ad esempio, la matrice A si trova sulla prima pagina d X e la sua decomposizione è data da U(:,:,1)*S(:,:,1)*V(:,:,1)'.

Approssimazioni SVD di rank ridotto

Per le grandi matrici sparse, non è sempre pratico utilizzare svd per calcolare tutti i valori singolari e i vettori singolari. Ad esempio, se fosse necessario conoscere solo alcuni dei valori singolari più grandi, calcolare tutti i valori singolari di una matrice sparsa 5000x5000 rappresenta un lavoro in più.

Nei casi in cui sia necessario solo un sottoinsieme dei valori singolari e dei vettori singolari, è preferibile utilizzare le funzioni svds e svdsketchrispetto a svd.

Funzione Utilizzo

svds Utilizzare svds per calcolare un'approssimazione rank-k della SVD. È possibile specificare se la sotto-serie di valori singolari deve essere quella più grande, quella più piccola, oppure quella più prossima a un numero specifico. svds generalmente calcola la miglior approssimazione possibile rank-k.

svdsketch Utilizzare svdsketch per calcolare un'approssimazione SVD parziale della matrice di input che soddisfa una tolleranza specifica. Mentre svds richiede la specifica del rank, svdsketch determina in modo adattivo il rank dello sketch della matrice in base alla tolleranza specificata. L'approssimazione rank-kutilizzata in ultima istanza da svdsketchsoddisfa la tolleranza ma, diversamente da svds, non è garantito che sia la migliore possibile.

Funzione	Utilizzo
`svds`	Utilizzare `svds` per calcolare un'approssimazione rank-k della SVD. È possibile specificare se la sotto-serie di valori singolari deve essere quella più grande, quella più piccola, oppure quella più prossima a un numero specifico. `svds` generalmente calcola la miglior approssimazione possibile rank-k.
`svdsketch`	Utilizzare `svdsketch` per calcolare un'approssimazione SVD parziale della matrice di input che soddisfa una tolleranza specifica. Mentre `svds` richiede la specifica del rank, `svdsketch` determina in modo adattivo il rank dello sketch della matrice in base alla tolleranza specificata. L'approssimazione rank-kutilizzata in ultima istanza da `svdsketch`soddisfa la tolleranza ma, diversamente da `svds`, non è garantito che sia la migliore possibile.

Ad esempio, si consideri una matrice sparsa casuale 1000x1000 con una densità di circa il 30%.

n = 1000;
A = sprand(n,n,0.3);

I sei valori singolari più grandi sono

S = svds(A)

S =

  130.2184
   16.4358
   16.4119
   16.3688
   16.3242
   16.2838

I sei valori singolari più piccoli, inoltre, sono

S = svds(A,6,"smallest")

S =

    0.0740
    0.0574
    0.0388
    0.0282
    0.0131
    0.0066

Per le matrici più piccole, che possono essere memorizzate come una matrice completa, full(A), l'uso di svd(full(A)) potrebbe risultare più veloce rispetto a svds o svdsketch. Per matrici molto grandi e sparse, tuttavia, diventa necessario usare svds o svdsketch.

Vedi anche

svd | svds | svdsketch | gsvd | pagesvd