Q
{\displaystyle Q}
è una matrice unitaria, cioè
Q
Q
h
=
Q
h
Q
=
1
{\displaystyle QQ^{h}=Q^{h}Q=1}
, allora
Q
−
1
=
Q
h
{\displaystyle Q^{-1}=Q^{h}}
, e
R
{\displaystyle R}
è una triangolare superiore.
Ci chiediamo per quale motivo serve un'alternativa alla fattorizzazione
P
A
=
L
U
{\displaystyle PA=LU}
.
Il vantaggio della fattorizzazione
Q
R
{\displaystyle QR}
è una maggior stabilità rispetto alla
P
A
=
L
U
{\displaystyle PA=LU}
. Infatti la
Q
R
{\displaystyle QR}
è una fattorizzazione con matrici elementari unitarie, che conservano la norma 2.
Lo svantaggio è il raddoppio del costo computazionale.
Non occorre la tecnica del pivot.
Ci sono situazioni più sofisticate in cui si usa una tecnica di pivot per colonne, per garantire che gli elementi diagonali della
R
{\displaystyle R}
abbiano una certa proprietà di ordinamento.
Questa fattorizzazione si usa solo in casi limite, ad esempio nel calcolo degli autovalori (spettro completo) e del rango, e della s.v.d.
Data una matrice
A
{\displaystyle A}
hermitiana, allora si può diagonalizzare con trasformazioni unitarie, scrivendo
Q
λ
Q
h
{\displaystyle Q\lambda Q^{h}}
.
La decomposizione in fattori singolari è l'estensione di quest'idea anche ai casi in cui la matrice non è simmetrica e non è quadrata.
La fattorizzazione
Q
R
{\displaystyle QR}
esiste sempre ma non è unica.
Con la fattorizzazione
Q
R
{\displaystyle QR}
A
{\displaystyle A}
viene fattorizzata nel prodotto di una matrice unitaria e di una triangolare superiore.
Q
{\displaystyle Q}
è unitaria quindi
Q
−
1
=
Q
h
{\displaystyle Q^{-1}=Q^{h}}
.
Nel sistema
A
x
=
b
{\displaystyle A{\boldsymbol {x}}={\boldsymbol {b}}}
, cioè
Q
R
x
=
b
{\displaystyle QR{\boldsymbol {x}}={\boldsymbol {b}}}
, poniamo
R
x
=
y
{\displaystyle R{\boldsymbol {x}}={\boldsymbol {y}}}
, risolvo prima
Q
y
=
b
⟶
y
=
Q
−
1
b
=
Q
h
b
prodotto matrice per vettore
{\displaystyle Q{\boldsymbol {y}}={\boldsymbol {b}}\,\longrightarrow {\boldsymbol {y}}=Q^{-1}{\boldsymbol {b}}=Q^{h}{\boldsymbol {b}}\quad {\hbox{prodotto matrice per vettore}}}
quindi si risolve
R
x
=
y
{\displaystyle R{\boldsymbol {x}}={\boldsymbol {y}}}
con
R
{\displaystyle R}
matrice triangolare.
Il costo di risoluzione è
2
/
3
n
3
{\displaystyle 2/3n^{3}}
.
Definizione
Una matrice elementare è una matrice della forma
E
=
I
−
σ
u
v
h
{\displaystyle E=I-\sigma {\boldsymbol {u}}{\boldsymbol {v}}^{h}}
, con
u
,
v
∈
C
n
{\displaystyle {\boldsymbol {u}},{\boldsymbol {v}}\in \mathbb {C} ^{n}}
,
σ
{\displaystyle \sigma }
parametro, e il prodotto
u
v
h
{\displaystyle {\boldsymbol {u}}{\boldsymbol {v}}^{h}}
è una diade.
La matrice della diade ha rango uguale a 1.
Nella fattorizzazione
Q
R
{\displaystyle QR}
consideriamo le matrici elementari di Hausolder (o matrici di riflessione). Queste matrici sono della forma:
P
=
I
−
β
v
v
h
,
β
∈
R
,
v
∈
C
n
{\displaystyle P=I-\beta {\boldsymbol {v}}{\boldsymbol {v}}^{h},\quad \beta \in \mathbb {R} ,\,{\boldsymbol {v}}\in \mathbb {C} ^{n}}
Le matrici di riflessione sono hermitiane, e sono unitarie se scegliamo
β
=
2
|
v
|
2
{\displaystyle \beta ={\frac {2}{\vert {\boldsymbol {v}}\vert _{2}}}}
.
Dim.
Hermitiane :
P
h
=
I
−
β
(
v
v
h
)
h
=
I
−
β
v
v
h
=
P
{\displaystyle P^{h}=I-\beta ({\boldsymbol {v}}{\boldsymbol {v}}^{h})^{h}=I-\beta {\boldsymbol {v}}{\boldsymbol {v}}^{h}=P}
Unitarie : sappiamo che
I
=
P
P
h
{\displaystyle I=PP^{h}}
, e per hermitianità
P
P
h
=
P
P
=
(
I
−
β
v
v
h
)
(
I
−
β
v
v
h
)
{\displaystyle PP^{h}=PP=(I-\beta {\boldsymbol {v}}{\boldsymbol {v}}^{h})(I-\beta {\boldsymbol {v}}{\boldsymbol {v}}^{h})}
=
I
−
β
v
v
h
−
β
v
v
h
+
β
2
v
v
h
v
v
h
{\displaystyle =I-\beta {\boldsymbol {v}}{\boldsymbol {v}}^{h}-\beta {\boldsymbol {v}}{\boldsymbol {v}}^{h}+\beta ^{2}{\boldsymbol {v}}{\boldsymbol {v}}^{h}{\boldsymbol {v}}{\boldsymbol {v}}^{h}}
=
I
−
2
β
v
v
h
+
β
2
v
(
|
v
|
)
2
v
h
{\displaystyle =I-2\beta {\boldsymbol {v}}{\boldsymbol {v}}^{h}+\beta ^{2}{\boldsymbol {v}}(\vert {\boldsymbol {v}}\vert )^{2}{\boldsymbol {v}}^{h}}
(infatti
v
h
v
=
(
|
v
|
2
)
2
{\displaystyle {\boldsymbol {v}}^{h}{\boldsymbol {v}}=(\vert {\boldsymbol {v}}\vert _{2})^{2}}
)
quindi
I
=
I
−
2
β
v
v
h
+
β
(
|
v
|
2
)
2
v
v
h
{\displaystyle I=I-2\beta {\boldsymbol {v}}{\boldsymbol {v}}^{h}+\beta (\vert {\boldsymbol {v}}\vert _{2})^{2}{\boldsymbol {v}}{\boldsymbol {v}}^{h}}
implica
β
v
v
h
∗
[
−
2
+
β
|
v
|
2
2
]
=
0
{\displaystyle \beta {\boldsymbol {v}}{\boldsymbol {v}}^{h}*[-2+\beta \vert {\boldsymbol {v}}\vert _{2}^{2}]=0}
cioè, se
β
≠
0
{\displaystyle \beta \neq 0}
, deve valere
−
2
+
β
|
v
|
2
2
=
0
{\displaystyle -2+\beta \vert {\boldsymbol {v}}\vert _{2}^{2}=0}
, cioè
β
=
2
(
|
v
|
2
)
2
{\displaystyle \beta ={\frac {2}{(\vert {\boldsymbol {v}}\vert _{2})^{2}}}}
cvd
Sia
P
=
P
(
v
)
{\displaystyle P=P({\boldsymbol {v}})}
(la matrice di Hausolder è univocamente determinata da
v
{\displaystyle {\boldsymbol {v}}}
), con
P
=
I
−
β
v
v
h
{\displaystyle P=I-\beta {\boldsymbol {v}}{\boldsymbol {v}}^{h}}
, e
β
=
2
(
|
v
|
2
)
2
{\displaystyle \beta ={\frac {2}{(\vert {\boldsymbol {v}}\vert _{2})^{2}}}}
. Allora per ogni
x
∈
v
⊥
{\displaystyle {\boldsymbol {x}}\in {\boldsymbol {v}}_{\perp }}
, la matrice di Hausolder si comporta come l'identità. Infatti:
x
(
I
−
β
v
v
h
)
=
x
−
β
v
v
h
x
=
x
{\displaystyle {\boldsymbol {x}}(I-\beta {\boldsymbol {v}}{\boldsymbol {v}}^{h})={\boldsymbol {x}}-\beta {\boldsymbol {v}}{\boldsymbol {v}}^{h}{\boldsymbol {x}}={\boldsymbol {x}}}
infatti
x
v
=
0
{\displaystyle {\boldsymbol {x}}{\boldsymbol {v}}=0}
per
x
∈
v
⊥
{\displaystyle {\boldsymbol {x}}\in {\boldsymbol {v}}_{\perp }}
.
P
{\displaystyle P}
applicata nel sottospazio
s
p
a
n
{
v
}
{\displaystyle \mathrm {span} \{{\boldsymbol {v}}\}}
è:
(
I
−
β
v
v
h
)
(
γ
v
)
=
γ
v
−
β
γ
v
v
h
v
{\displaystyle (I-\beta {\boldsymbol {v}}{\boldsymbol {v}}^{h})(\gamma {\boldsymbol {v}})=\gamma {\boldsymbol {v}}-\beta \gamma {\boldsymbol {v}}{\boldsymbol {v}}^{h}{\boldsymbol {v}}}
=
γ
v
−
β
γ
|
v
|
2
v
{\displaystyle =\gamma {\boldsymbol {v}}-\beta \gamma \vert {\boldsymbol {v}}\vert ^{2}{\boldsymbol {v}}}
e per definizione di
β
{\displaystyle \beta }
:
=
γ
v
−
2
(
|
v
|
2
)
2
γ
|
v
|
2
v
{\displaystyle =\gamma {\boldsymbol {v}}-{\frac {2}{(\vert {\boldsymbol {v}}\vert _{2})^{2}}}\gamma \vert {\boldsymbol {v}}\vert ^{2}{\boldsymbol {v}}}
=
γ
v
−
2
γ
v
=
−
γ
v
{\displaystyle =\gamma {\boldsymbol {v}}-2\gamma {\boldsymbol {v}}=-\gamma {\boldsymbol {v}}}
cioè si ha una riflessione.
Per ogni
x
∈
C
n
∖
0
{\displaystyle {\boldsymbol {x}}\in \mathbb {C} ^{n}\smallsetminus 0}
esiste
P
∈
C
n
×
n
{\displaystyle P\in \mathbb {C} _{n\times n}}
matrice elementare di Hausolder tale che
P
x
=
α
e
1
{\displaystyle P{\boldsymbol {x}}=\alpha {\boldsymbol {e}}_{1}}
.
Dobbiamo determinare il parametro
α
∈
C
{\displaystyle \alpha \in \mathbb {C} }
e la matrice elementare di Hausolder
P
{\displaystyle P}
.
Siccome
α
∈
C
{\displaystyle \alpha \in \mathbb {C} }
,
α
=
|
α
|
∗
e
i
θ
{\displaystyle \alpha =|\alpha |*e^{i\theta }}
(relazione
∘
{\displaystyle \circ }
). Deve essere
P
x
=
α
e
1
{\displaystyle P{\boldsymbol {x}}=\alpha {\boldsymbol {e}}_{1}}
, allora
|
P
x
|
=
|
α
|
|
e
1
|
2
=
|
α
|
{\displaystyle \vert P{\boldsymbol {x}}\vert =|\alpha |\vert {\boldsymbol {e}}_{1}\vert _{2}=|\alpha |}
Inoltre
P
{\displaystyle P}
è un'isometria quindi conserva la norma 2, cioè:
(
|
P
x
|
2
)
2
=
P
2
x
h
x
=
x
h
x
=
(
|
x
|
2
)
2
{\displaystyle (\vert P{\boldsymbol {x}}\vert _{2})^{2}=P^{2}{\boldsymbol {x}}^{h}{\boldsymbol {x}}={\boldsymbol {x}}^{h}{\boldsymbol {x}}=(\vert {\boldsymbol {x}}\vert _{2})^{2}}
allora, eguagliando le due espressioni trovate per
(
|
P
x
|
2
)
2
{\displaystyle (\vert P{\boldsymbol {x}}\vert _{2})^{2}}
, ottengo:
|
α
|
=
|
x
|
2
{\displaystyle |\alpha |=\vert {\boldsymbol {x}}\vert _{2}}
Per determinare l'angolo
θ
{\displaystyle \theta }
, dev'essere
x
h
P
x
=
x
h
α
e
1
=
α
x
h
e
1
=
α
x
¯
1
è uno scalare
{\displaystyle {\boldsymbol {x}}^{h}P{\boldsymbol {x}}={\boldsymbol {x}}^{h}\alpha {\boldsymbol {e}}_{1}=\alpha {\boldsymbol {x}}^{h}{\boldsymbol {e}}_{1}=\alpha {\bar {x}}_{1}\quad {\hbox{ è uno scalare}}}
Per hermitianità di
P
{\displaystyle P}
:
(
¯
x
h
P
x
)
=
(
x
h
P
x
)
h
=
x
h
P
h
x
=
x
h
P
x
{\displaystyle {\bar {(}}{\boldsymbol {x}}^{h}P{\boldsymbol {x}})=({\boldsymbol {x}}^{h}P{\boldsymbol {x}})^{h}={\boldsymbol {x}}^{h}P^{h}{\boldsymbol {x}}={\boldsymbol {x}}^{h}P{\boldsymbol {x}}}
allora siccome lo scalare è uguale al suo coniugato, è reale, quindi
α
x
¯
1
∈
R
{\displaystyle \alpha {\bar {x}}_{1}\in \mathbb {R} }
.
x
¯
1
∈
C
{\displaystyle {\bar {x}}_{1}\in \mathbb {C} }
, allora possiamo scriverlo come
|
x
1
|
e
i
ϕ
{\displaystyle |x_{1}|e^{i\phi }}
(relazione
∘
∘
{\displaystyle \circ \circ }
).
α
x
¯
1
=
|
α
|
e
i
θ
|
x
1
|
e
i
ϕ
{\displaystyle \alpha {\bar {x}}_{1}=|\alpha |e^{i\theta }|x_{1}|e^{i\phi }}
e sostituendo l'espressione di
α
{\displaystyle \alpha }
:
=
|
x
|
2
|
x
1
|
e
i
(
θ
−
ϕ
)
{\displaystyle =\vert {\boldsymbol {x}}\vert _{2}|x_{1}|e^{i(\theta -\phi )}}
α
x
¯
1
∈
R
{\displaystyle \alpha {\bar {x}}_{1}\in \mathbb {R} }
, allora
θ
−
ϕ
=
k
π
{\displaystyle \theta -\phi =k\pi }
ovvero
θ
=
ϕ
+
δ
π
{\displaystyle \theta =\phi +\delta \pi }
con
δ
=
0
,
1
{\displaystyle \delta =0,1}
.
In conclusione:
α
=
|
α
|
e
i
θ
=
|
x
|
2
∗
e
i
(
ϕ
+
δ
π
)
{\displaystyle \alpha =|\alpha |e^{i\theta }=\vert {\boldsymbol {x}}\vert _{2}*e^{i(\phi +\delta \pi )}}
e
i
δ
π
=
±
1
⟶
α
=
±
|
x
|
2
e
i
ϕ
{\displaystyle e^{i\delta \pi }=\pm 1\,\longrightarrow \,\alpha =\pm \vert {\boldsymbol {x}}\vert _{2}e^{i\phi }}
e per la
∘
∘
{\displaystyle \circ \circ }
,
=
±
|
x
|
2
∗
x
1
|
x
1
|
{\displaystyle =\pm \vert {\boldsymbol {x}}\vert _{2}*{\frac {x_{1}}{|x_{1}|}}}
Se
x
1
=
0
{\displaystyle x_{1}=0}
, allora poniamo
α
=
−
|
x
|
2
{\displaystyle \alpha =-\vert {\boldsymbol {x}}\vert _{2}}
.
Vale che
P
x
=
(
I
−
β
v
v
h
)
(
x
)
=
x
−
β
v
v
h
x
,
{\displaystyle P{\boldsymbol {x}}=(I-\beta {\boldsymbol {v}}{\boldsymbol {v}}^{h})({\boldsymbol {x}})={\boldsymbol {x}}-\beta {\boldsymbol {v}}{\boldsymbol {v}}^{h}{\boldsymbol {x}},}
e per definizione di
β
{\displaystyle \beta }
otteniamo:
P
x
=
x
−
2
(
|
v
|
2
)
2
v
h
x
v
{\displaystyle P{\boldsymbol {x}}={\boldsymbol {x}}-{\frac {2}{(\vert {\boldsymbol {v}}\vert _{2})^{2}}}{\boldsymbol {v}}^{h}{\boldsymbol {x}}{\boldsymbol {v}}}
Chiamiamo
β
^
{\displaystyle {\hat {\beta }}}
il coefficiente che moltiplica
v
{\displaystyle {\boldsymbol {v}}}
, e mettiamo in evidenza
β
^
v
{\displaystyle {\hat {\beta }}{\boldsymbol {v}}}
:
β
^
v
=
x
−
P
x
=
x
−
α
e
1
{\displaystyle {\hat {\beta }}{\boldsymbol {v}}={\boldsymbol {x}}-P{\boldsymbol {x}}={\boldsymbol {x}}-\alpha {\boldsymbol {e}}_{1}}
A patto di escludere che
β
^
=
0
{\displaystyle {\hat {\beta }}=0}
, è lecito porre
β
^
=
1
{\displaystyle {\hat {\beta }}=1}
, in quanto
P
(
γ
v
)
=
I
−
2
|
γ
v
|
2
∗
(
γ
v
)
∗
(
γ
v
h
)
{\displaystyle P(\gamma {\boldsymbol {v}})=I-{\frac {2}{\vert \gamma {\boldsymbol {v}}\vert ^{2}}}*(\gamma {\boldsymbol {v}})*(\gamma {\boldsymbol {v}}^{h})}
=
I
−
2
|
γ
|
2
|
v
|
2
|
γ
|
2
∗
v
v
h
{\displaystyle =I-{\frac {2}{|\gamma |^{2}\vert {\boldsymbol {v}}\vert _{2}}}|\gamma |^{2}*{\boldsymbol {v}}{\boldsymbol {v}}^{h}}
=
I
−
2
|
v
|
2
∗
v
v
h
=
P
(
v
)
{\displaystyle =I-{\frac {2}{\vert {\boldsymbol {v}}\vert _{2}}}*{\boldsymbol {v}}{\boldsymbol {v}}_{h}=P({\boldsymbol {v}})}
allora per il fatto che
P
(
γ
v
)
=
P
(
v
)
{\displaystyle P(\gamma {\boldsymbol {v}})=P({\boldsymbol {v}})}
, possiamo porre
β
^
=
1
{\displaystyle {\hat {\beta }}=1}
.
Tornando alla formula:
β
^
v
=
x
−
α
e
1
{\displaystyle {\hat {\beta }}{\boldsymbol {v}}={\boldsymbol {x}}-\alpha {\boldsymbol {e}}_{1}}
e ponendo
β
^
=
1
{\displaystyle {\hat {\beta }}=1}
, otteniamo:
v
=
x
−
α
e
1
{\displaystyle {\boldsymbol {v}}={\boldsymbol {x}}-\alpha {\boldsymbol {e}}_{1}}
Nota bene: possiamo escludere il caso
β
^
=
0
{\displaystyle {\hat {\beta }}=0}
, infatti, per definizione,
β
^
=
0
{\displaystyle {\hat {\beta }}=0}
se e solo se
x
∈
v
⊥
{\displaystyle {\boldsymbol {x}}\in {\boldsymbol {v}}_{\perp }}
, e possiamo scegliere il vettore
x
{\displaystyle {\boldsymbol {x}}}
al di fuori del complemento ortogonale.
Calcoliamo:
v
h
v
=
v
h
(
x
−
α
e
1
)
=
(
|
x
|
)
2
+
|
x
1
|
|
x
|
2
{\displaystyle {\boldsymbol {v}}^{h}{\boldsymbol {v}}={\boldsymbol {v}}^{h}({\boldsymbol {x}}-\alpha {\boldsymbol {e}}_{1})=(\vert {\boldsymbol {x}}\vert )^{2}+|x_{1}|\vert {\boldsymbol {x}}\vert _{2}}
e la quantità è positiva se
x
≠
0
{\displaystyle {\boldsymbol {x}}\neq 0}
come nell'ipotesi.
cvd
Riassumendo, poniamo
α
=
{
±
|
x
|
2
x
1
|
x
1
|
se
x
1
≠
0
−
|
x
|
2
se
x
1
=
0
{\displaystyle \alpha ={\begin{cases}\pm \vert {\boldsymbol {x}}\vert _{2}{\frac {x_{1}}{|x_{1}|}}&{\hbox{se}}x_{1}\neq 0\\-\vert {\boldsymbol {x}}\vert _{2}&{\hbox{se}}x_{1}=0\end{cases}}}
v
=
x
−
α
e
1
{\displaystyle {\boldsymbol {v}}={\boldsymbol {x}}-\alpha {\boldsymbol {e}}_{1}}
e
1
{\displaystyle {\boldsymbol {e}}_{1}}
modifica solo la prima componente del vettore
v
{\displaystyle {\boldsymbol {v}}}
. Quindi:
v
1
=
x
1
−
α
,
=
x
1
−
(
−
|
x
|
2
x
1
|
x
1
|
)
=
{\displaystyle v_{1}=x_{1}-\alpha ,=x_{1}-(-\vert {\boldsymbol {x}}\vert _{2}{\frac {x_{1}}{|x_{1}|}})=}
=
x
1
+
|
x
|
2
x
1
|
x
1
|
)
=
{\displaystyle =x_{1}+\vert {\boldsymbol {x}}\vert _{2}{\frac {x_{1}}{|x_{1}|}})=}
=
x
1
∗
(
|
x
1
|
+
|
x
|
2
|
x
1
|
)
{\displaystyle =x_{1}*({\frac {|x_{1}|+\vert {\boldsymbol {x}}\vert _{2}}{|x_{1}|}})}
e non c'è possibilità di errori di cancellazione.
Se
x
∈
R
n
{\displaystyle {\boldsymbol {x}}\in \mathbb {R} ^{n}}
,
α
∈
R
,
v
∈
R
n
{\displaystyle \alpha \in \mathbb {R} ,{\boldsymbol {v}}\in \mathbb {R} ^{n}}
per costruzione, allora
P
{\displaystyle P}
è una simmetrica ortogonale.
Il calcolo di
v
{\displaystyle {\boldsymbol {v}}}
costa
n
{\displaystyle n}
operazioni per il calcolo della norma di
x
{\displaystyle {\boldsymbol {x}}}
, 2 operazioni moltiplicative, un'estrazione di radice, cioè un numero dell'ordine di
n
{\displaystyle n}
operazioni.
Utilizziamo il teorema dimostrato per determinare la fattorizzazione
Q
R
{\displaystyle QR}
.
Supponiamo di essere al passo
k
{\displaystyle k}
-esimo della fattorizzazione
Q
R
{\displaystyle QR}
. Allora chiamiamo
a
k
=
(
A
k
k
,
…
,
A
n
k
)
{\displaystyle {\boldsymbol {a}}_{k}=(A_{kk},\dots ,A_{nk})}
la
k
{\displaystyle k}
-esima sottocolonna della matrice
A
(
k
)
{\displaystyle A^{(k)}}
, che dobbiamo azzerare (corrisponde al vettore
x
{\displaystyle {\boldsymbol {x}}}
del teorema). Scegliamo il vettore
v
k
=
v
{\displaystyle {\boldsymbol {v}}_{k}={\boldsymbol {v}}}
che ha 0 nelle prime
k
−
1
{\displaystyle k-1}
componenti, e le componenti successive sono:
v
k
=
γ
k
(
|
a
k
|
2
+
|
A
k
k
(
k
)
|
)
=
s
g
n
a
k
k
∗
(
|
σ
k
|
)
{\displaystyle {\boldsymbol {v}}_{k}=\gamma _{k}(\vert {\boldsymbol {a}}_{k}\vert _{2}+|A_{kk}^{(k)}|)=\mathrm {sgn} a_{kk}*(\vert \sigma _{k}\vert )}
v
k
+
1
=
a
k
+
1
(
k
)
{\displaystyle {\boldsymbol {v}}_{k+1}=a_{k+1}^{(k)}}
v
k
+
2
=
a
k
+
2
(
k
)
{\displaystyle {\boldsymbol {v}}_{k+2}=a_{k+2}^{(k)}}
…
{\displaystyle \dots }
v
n
=
a
n
(
k
)
{\displaystyle {\boldsymbol {v}}_{n}=a_{n}^{(k)}}
dove definiamo
γ
k
=
{
A
k
k
|
A
k
k
|
⟺
A
k
k
≠
0
1
se
A
k
k
=
0
{\displaystyle \gamma _{k}={\begin{cases}{\frac {A_{kk}}{|A_{kk}|}}&\iff A_{kk}\neq 0\\1&{\hbox{se}}A_{kk}=0\end{cases}}}
(per la definizione di
v
k
{\displaystyle {\boldsymbol {v}}_{k}}
sfruttiamo l'espressione ricavata nel teorema, cioè
v
k
=
x
−
α
e
1
{\displaystyle {\boldsymbol {v}}_{k}={\boldsymbol {x}}-\alpha {\boldsymbol {e}}_{1}}
, e sfruttiamo anche l'espressione ricavata per
v
1
{\displaystyle v_{1}}
)
β
k
=
2
(
|
v
k
|
2
)
2
{\displaystyle \beta _{k}={\frac {2}{(\vert {\boldsymbol {v}}_{k}\vert _{2})^{2}}}}
Allora calcoliamo:
|
v
k
|
2
=
|
γ
k
|
2
(
|
a
k
|
2
+
|
A
k
k
|
)
2
+
|
a
k
|
2
2
−
|
A
k
k
|
2
{\displaystyle \vert {\boldsymbol {v}}_{k}\vert ^{2}=|\gamma _{k}|^{2}(\vert {\boldsymbol {a}}_{k}\vert _{2}+|A_{kk}|)^{2}+\vert {\boldsymbol {a}}_{k}\vert _{2}^{2}-|A_{kk}|^{2}}
=
(
|
a
k
|
2
)
2
+
2
|
A
k
k
|
|
a
k
|
2
+
(
|
A
k
k
|
)
2
+
|
a
k
|
2
2
−
|
A
k
k
|
2
{\displaystyle =(\vert {\boldsymbol {a}}_{k}\vert _{2})^{2}+2|A_{kk}|\vert {\boldsymbol {a}}_{k}\vert _{2}+(|A_{kk}|)^{2}+\vert {\boldsymbol {a}}_{k}\vert _{2}^{2}-|A_{kk}|^{2}}
=
2
(
|
a
k
|
2
)
2
+
2
|
A
k
k
|
|
a
k
|
2
{\displaystyle =2(\vert {\boldsymbol {a}}_{k}\vert _{2})^{2}+2|A_{kk}|\vert {\boldsymbol {a}}_{k}\vert _{2}}
e raccogliendo
2
|
a
k
|
2
{\displaystyle 2\vert {\boldsymbol {a}}_{k}\vert _{2}}
otteniamo:
=
2
|
a
k
|
2
∗
[
|
a
k
|
2
+
|
A
k
k
|
]
{\displaystyle =2\vert {\boldsymbol {a}}_{k}\vert _{2}*[\vert {\boldsymbol {a}}_{k}\vert _{2}+|A_{kk}|]}
cioè
β
k
=
2
2
(
|
v
k
|
2
)
2
{\displaystyle \beta _{k}={\frac {2}{2(\vert {\boldsymbol {v}}_{k}\vert _{2})^{2}}}}
=
2
2
|
a
k
|
2
∗
[
|
a
k
|
+
A
k
k
]
{\displaystyle ={\frac {2}{2\vert {\boldsymbol {a}}_{k}\vert _{2}*[\vert {\boldsymbol {a}}_{k}\vert +A_{kk}]}}}
β
k
=
1
|
a
k
|
2
∗
[
|
a
k
|
+
A
k
k
]
{\displaystyle \beta _{k}={\frac {1}{\vert {\boldsymbol {a}}_{k}\vert _{2}*[\vert {\boldsymbol {a}}_{k}\vert +A_{kk}]}}}
Allora abbiamo determinato univocamente la matrice
P
k
=
I
−
β
v
v
h
{\displaystyle P_{k}=I-\beta {\boldsymbol {v}}{\boldsymbol {v}}^{h}}
Se il vettore
a
k
{\displaystyle {\boldsymbol {a}}_{k}}
è tutto nullo, consideriamo
γ
k
=
1
{\displaystyle \gamma _{k}=1}
e
v
k
=
0
{\displaystyle {\boldsymbol {v}}_{k}=0}
, in modo che a quel passo la matrice di Hausolder sia l'identità.
Applicando la matrice
P
{\displaystyle P}
trovata al vettore
x
{\displaystyle {\boldsymbol {x}}}
cioè alla colonna
a
k
{\displaystyle {\boldsymbol {a}}_{k}}
, tutte le entrate di
a
k
{\displaystyle {\boldsymbol {a}}_{k}}
da
k
+
1
{\displaystyle k+1}
a
n
{\displaystyle n}
vengono azzerate, mentre si ha
v
k
=
A
k
k
=
−
|
x
|
2
{\displaystyle v_{k}=A_{kk}=-\vert {\boldsymbol {x}}\vert _{2}}
(infatti
P
x
=
α
e
1
{\displaystyle P{\boldsymbol {x}}=\alpha {\boldsymbol {e}}_{1}}
). Applichiamo poi la matrice di Hausolder alle colonne rimanenti di
A
{\displaystyle A}
.
Riepilogando, data la matrice
A
{\displaystyle A}
di partenza, con il primo passo
P
1
A
{\displaystyle P_{1}A}
azzero la prima colonna e modifico
A
11
{\displaystyle A_{11}}
, e poi ripeto il procedimento, per
n
−
1
{\displaystyle n-1}
passi, e otteniamo la matrice triangolare superiore
R
{\displaystyle R}
.
Il prodotto di matrici unitarie è unitario, e sappiamo che
R
=
P
n
−
1
…
P
2
P
1
A
{\displaystyle R=P_{n-1}\dots P_{2}P_{1}A}
e per hermitianità:
A
=
P
1
P
2
…
P
n
−
1
R
{\displaystyle A=P_{1}P_{2}\dots P_{n-1}R}
quindi
Q
=
P
1
P
2
…
P
n
−
1
{\displaystyle Q=P_{1}P_{2}\dots P_{n-1}}
Per ogni
A
{\displaystyle A}
esiste la fattorizzazione
A
=
Q
R
{\displaystyle A=QR}
, ma non è unica.
Dim. Consideriamo una matrice di fase
S
{\displaystyle S}
, cioè una matrice diagonale con elementi
θ
1
,
…
,
θ
n
{\displaystyle \theta _{1},\dots ,\theta _{n}}
con
|
θ
i
|
=
1
∀
i
{\displaystyle |\theta _{i}|=1\forall i}
, ed è una matrice unitaria. Allora
S
S
h
=
I
{\displaystyle SS_{h}=I}
, cioè
A
=
Q
R
=
Q
S
S
h
R
{\displaystyle A=QR=QSS_{h}R}
, e definiamo
Q
¯
,
R
¯
{\displaystyle {\bar {Q}},{\bar {R}}}
dove
Q
¯
=
Q
S
{\displaystyle {\bar {Q}}=QS}
e
R
¯
=
S
h
R
{\displaystyle {\bar {R}}=S_{h}R}
, dove
Q
¯
{\displaystyle {\bar {Q}}}
è ancora unitaria e
R
¯
{\displaystyle {\bar {R}}}
è triangolare superiore. Allora abbiamo comunque una fattorizzazione
Q
R
{\displaystyle QR}
della matrice.
cvd
Questa proprietà può essere sfruttata per scegliere la
R
¯
{\displaystyle {\bar {R}}}
in modo che gli elementi sulla diagonale siano ordinati in una certa maniera.
Per risolvere il sistema lineare, non è necessario trovare esplicitamente la matrice
Q
{\displaystyle Q}
, e nemmeno le
P
h
{\displaystyle P_{h}}
.
Passo 1 : poniamo
A
1
=
A
,
b
1
=
b
{\displaystyle A_{1}=A,{\boldsymbol {b}}_{1}={\boldsymbol {b}}}
, e chiamiamo
T
1
{\displaystyle T_{1}}
la matrice che accosta
b
1
{\displaystyle {\boldsymbol {b}}_{1}}
.
Passo 2 : chiamiamo
T
2
=
P
1
T
1
{\displaystyle T_{2}=P_{1}T_{1}}
, con
P
1
=
I
−
β
1
v
1
v
1
h
{\displaystyle P_{1}=I-\beta _{1}{\boldsymbol {v}}_{1}{\boldsymbol {v}}_{1}^{h}}
Chiamiamo poi
y
1
h
{\displaystyle {\boldsymbol {y}}_{1}^{h}}
il vettore riga
v
1
h
T
1
{\displaystyle {\boldsymbol {v}}_{1}^{h}T_{1}}
, in modo da scrivere:
T
2
=
P
1
T
1
=
T
1
−
β
1
v
1
y
1
h
{\displaystyle T_{2}=P_{1}T_{1}=T_{1}-\beta _{1}{\boldsymbol {v}}_{1}{\boldsymbol {y}}_{1}^{h}}
Passo
k
+
1
{\displaystyle k+1}
:
T
k
+
1
=
P
k
T
k
=
T
k
−
β
k
v
k
y
k
h
{\displaystyle T_{k+1}=P_{k}T_{k}=T_{k}-\beta _{k}{\boldsymbol {v}}_{k}{\boldsymbol {y}}_{k}^{h}}
con
y
k
h
=
v
k
h
T
k
{\displaystyle {\boldsymbol {y}}_{k}^{h}={\boldsymbol {v}}_{k}^{h}T_{k}}
Infine
Passo
n
{\displaystyle n}
:
T
n
=
R
⟶
y
k
h
=
v
k
R
{\displaystyle T_{n}=R\,\longrightarrow {\boldsymbol {y}}_{k}^{h}={\boldsymbol {v}}_{k}R}
il vettore
y
{\displaystyle {\boldsymbol {y}}}
è stato calcolato senza calcolare esplicitamente le matrici
Q
{\displaystyle Q}
, ed è tale che
q
y
=
b
{\displaystyle q{\boldsymbol {y}}={\boldsymbol {b}}}
.
Per analizzare il costo del
k
{\displaystyle k}
-esimo passo, tengo conto che
v
k
{\displaystyle {\boldsymbol {v}}_{k}}
ha
k
−
1
{\displaystyle k-1}
componenti nulle e
n
−
k
+
1
{\displaystyle n-k+1}
componenti diverse da 0. Per calcolare
v
k
{\displaystyle {\boldsymbol {v}}_{k}}
, le operazioni sono quindi
n
−
k
+
1
+
2
=
n
−
k
+
3
{\displaystyle n-k+1+2=n-k+3}
.
…
…
…
{\displaystyle \dots \dots \dots }
Definizione
Si definisce norma di Frobenius' Testo in grassetto' di
A
{\displaystyle A}
la quantità
∑
i
,
j
=
1
n
|
A
i
j
|
2
{\displaystyle {\sqrt {\sum _{i,j=1}^{n}|A_{ij}|^{2}}}}
=
t
r
(
A
h
∗
A
)
{\displaystyle ={\sqrt {\mathrm {tr} (A_{h}*A)}}}
Questa non è una norma matriciale indotta, infatti
|
I
|
f
=
n
{\displaystyle \vert I\vert _{f}={\sqrt {n}}}
ed è diversa da 1, mentre con le norme matriciali indotte l'identità a norma 1.
Sia
x
¯
{\displaystyle {\bar {\boldsymbol {x}}}}
la soluzione effettivamente calcolata con l'applicazione delle
P
¯
k
{\displaystyle {\bar {P}}_{k}}
effettive alla matrice
A
{\displaystyle A}
per calcolare la
R
¯
{\displaystyle {\bar {R}}}
. Sia
y
¯
{\displaystyle {\bar {\boldsymbol {y}}}}
il vettore che compare nella risoluzione backward del sistema
R
x
=
y
{\displaystyle R{\boldsymbol {x}}={\boldsymbol {y}}}
. Allora
x
¯
{\displaystyle {\bar {\boldsymbol {x}}}}
è soluzione esatta di un sistema
(
A
+
δ
A
)
x
¯
=
b
+
δ
b
{\displaystyle (A+\delta _{A}){\bar {\boldsymbol {x}}}={\boldsymbol {b}}+\delta _{b}}
dove
|
δ
A
|
f
<
u
γ
n
2
|
A
|
f
+
n
|
R
¯
|
f
+
o
(
u
2
)
{\displaystyle \vert \delta _{A}\vert _{f}<u\gamma n^{2}\vert A\vert _{f}+n\vert {\bar {R}}\vert _{f}+o(u^{2})}
|
δ
b
|
2
<
u
γ
n
2
|
b
|
2
+
O
(
u
2
)
{\displaystyle \vert \delta _{b}\vert _{2}<u\gamma n^{2}\vert {\boldsymbol {b}}\vert _{2}+O(u^{2})}
L'entità della perturbazione dipende dai dati iniziali, dalle dimensioni della matrice, e dalla matrice di fattorizzazione. Non dipende dalla norma di Frobenius di
Q
{\displaystyle Q}
che è 1.
Consideriamo le matrici
A
(
k
)
{\displaystyle A^{(k)}}
al
k
{\displaystyle k}
-esimo passo di fattorizzazione, e calcolo il prodotto
A
k
+
1
h
A
k
+
1
=
A
k
h
P
k
h
P
k
A
k
{\displaystyle A_{k+1}^{h}A_{k+1}=A_{k}^{h}P_{k}^{h}P_{k}A_{k}}
siccome
P
k
h
P
k
=
I
{\displaystyle P_{k}^{h}P_{k}=I}
allora
A
k
h
A
k
{\displaystyle A_{k}^{h}A_{k}}
non varia con i passi di fattorizzazione.
Allora
t
r
(
A
k
+
1
h
A
k
+
1
)
=
t
r
(
A
k
h
A
k
)
{\displaystyle \mathrm {tr} (A_{k+1}^{h}A_{k+1})=\mathrm {tr} (A_{k}^{h}A_{k})}
⟶
|
A
k
+
1
|
f
2
=
|
A
k
|
f
2
{\displaystyle \longrightarrow \vert A_{k+1}\vert _{f}^{2}=\vert A_{k}\vert _{f}^{2}}
allora
|
A
k
|
f
=
|
A
|
f
{\displaystyle \vert A_{k}\vert _{f}=\vert A\vert _{f}}
inoltre
|
R
|
f
=
|
A
n
|
f
=
|
A
|
f
{\displaystyle \vert R\vert _{f}=\vert A_{n}\vert _{f}=\vert A\vert _{f}}
allora non c'è possibilità di esplosione dei coefficienti.
La
Q
R
{\displaystyle QR}
è stabile e può essere fatta senza pivot.