RDODT and RUPDT Calls
downdate and update QR and Cholesky decompositions
- CALL RDODT( def, rup, bup, sup,
,
,
,
,
ssq>);
- CALL RUPDT( rup, bup, sup,
,
,
,
,
ssq>);
The RDODT and RUPDT subroutines return the values:
- def
- is only used for downdating, and it specifies
whether the downdating of matrix
by using
the
rows in argument
has been successful.
The result def=2 means that the downdating
of
by at least one row of
leads to a
singular matrix and cannot be completed successfully
(since the result of downdating is not unique).
In that case, the results rup, bup,
and sup contain missing values only.
The result def=1 means that the residual sum of
squares, ssq, could not be downdated successfully
and the result sup contains missing values only.
The result def=0 means that the downdating
of
by
was completed successfully.
- rup
- is the
upper triangular matrix
that has
been updated or downdated by using the
rows in
.
- bup
- is the
matrix
of right-hand sides that has
been updated or downdated by using the
rows in argument
.
If the argument
is not specified, bup is not computed.
- sup
- is a
vector of square roots of residual sum of squares that
is updated or downdated by using the
rows of argument
.
If ssq is not specified, sup is not computed.
The inputs to the RDODT and
RUPDT subroutines are as follows:
![r](images/langref_langrefeq50.gif)
- specifies an
upper triangular matrix
to be updated or downdated by the
rows in
.
Only the upper triangle of
is used; the
lower triangle can contain any information.
![z](images/langref_langrefeq671.gif)
- specifies a
matrix
used
rowwise to update or downdate the matrix
.
![b](images/langref_langrefeq43.gif)
- specifies an optional
matrix
of right-hand sides that have to be
updated or downdated simultaneously with
.
If
is specified, the argument
must also be specified.
![y](images/langref_langrefeq145.gif)
- specifies an optional
matrix
used rowwise
to update or downdate the right-hand side matrix
.
If
is specified, the argument
must also be specified.
- ssq
- is an optional
vector that, if
is specified, specifies
the square root of the error sum of squares that should be
updated or downdated simultaneously with
and
.
The upper triangular matrix
![{{r}}](images/langref_langrefeq40.gif)
of the QR
decomposition of an
![m x n](images/langref_langrefeq34.gif)
matrix
![{a}](images/langref_langrefeq35.gif)
,
![{a}= {q}{{r}}, { where } {q}^' {q}= {q}{q}^' = {i}_m](images/langref_langrefeq1076.gif)
is recomputed efficiently in two cases:
- update: An
vector
is added to matrix
. - downdate: An
vector
is deleted from matrix
.
Computing the whole QR decomposition of matrix
![{a}](images/langref_langrefeq35.gif)
by
Householder transformations requires
![4mn^2 - 4n^3/3](images/langref_langrefeq1078.gif)
floating-point operations, whereas updating or downdating
the QR decomposition (by Givens rotations) of one row
vector
![{z}](images/langref_langrefeq1077.gif)
requires only
![2n^2](images/langref_langrefeq441.gif)
floating-point operations.
If the QR decomposition is used to solve
the full-rank linear least squares problem
![\min_{{x}} \vert {a}{x}- {b}\vert^2 = {ssq}](images/langref_langrefeq1079.gif)
by solving the nonsingular upper triangular system
![{x}= {{r}}^{-1} {q}^' {b}](images/langref_langrefeq1080.gif)
then the
RUPDT and RDODT subroutines can be used to update or
downdate the
![p](images/langref_langrefeq12.gif)
-transformed right-hand sides
![{q}^' {b}](images/langref_langrefeq149.gif)
and the residual sum-of-squares
![p](images/langref_langrefeq12.gif)
vector
ssq
provided that for each
![n](images/langref_langrefeq53.gif)
vector
![{z}](images/langref_langrefeq1077.gif)
added to or deleted
from
![{a}](images/langref_langrefeq35.gif)
there is also a
![p](images/langref_langrefeq12.gif)
vector
![{y}](images/langref_langrefeq758.gif)
added to or
deleted from the
![m x p](images/langref_langrefeq44.gif)
right-hand-side matrix
![{b}](images/langref_langrefeq45.gif)
.
If the arguments
![z](images/langref_langrefeq671.gif)
and
![y](images/langref_langrefeq145.gif)
of the subroutines
RUPDT and
RDODT contain
![q \gt 1](images/langref_langrefeq1081.gif)
row vectors for which
![{{r}}](images/langref_langrefeq40.gif)
(and
![{q}^' {b}](images/langref_langrefeq149.gif)
, and eventually
ssq) is
to be updated or downdated, the process is performed
stepwise by processing the rows
![{z}_k](images/langref_langrefeq1082.gif)
(and
![{y}_k](images/langref_langrefeq1083.gif)
),
![k = 1, ... ,q](images/langref_langrefeq1084.gif)
, in the order in which they are stored.
The QR decomposition of an
![m x n](images/langref_langrefeq34.gif)
matrix
![{a}](images/langref_langrefeq35.gif)
,
![m \geq n](images/langref_langrefeq36.gif)
, rank
![({a}) = n](images/langref_langrefeq1085.gif)
,
![{a}= {q}{{r}}, { where } {q}^' {q}= {q}{q}^' = {i}_m](images/langref_langrefeq1086.gif)
corresponds to the Cholesky factorization
![{c}= {{r}}^' {{r}}, { where } {c}= {a}^' {a}](images/langref_langrefeq1087.gif)
of the positive definite
![n x n](images/langref_langrefeq39.gif)
crossproduct matrix
![{c}= {a}^' {a}](images/langref_langrefeq1088.gif)
.
In the case where
![m \geq n](images/langref_langrefeq36.gif)
and rank
![({a}) = n](images/langref_langrefeq1085.gif)
, the upper
triangular matrix
![{{r}}](images/langref_langrefeq40.gif)
computed by the QR decomposition
(with positive diagonal elements) is the same as the one
computed by Cholesky factorization except for numerical error,
![{a}^' {a}= ({q}{{r}})^' ({q}{{r}}) = {{r}}^' {{r}}](images/langref_langrefeq1089.gif)
Adding a row vector
![{z}](images/langref_langrefeq1077.gif)
to matrix
![{a}](images/langref_langrefeq35.gif)
corresponds to
the rank-1 modification of the crossproduct matrix
![{\widetilde{c}} = {c}+ {z}^' {z}, { where } {\widetilde{c}} = {\widetilde{a}^' \widetilde{a}}](images/langref_langrefeq1091.gif)
and the
![(m+1) x n](images/langref_langrefeq1092.gif)
matrix
![{\widetilde{a}}](images/langref_langrefeq1093.gif)
contains all rows of
![{a}](images/langref_langrefeq35.gif)
with the row
![{z}](images/langref_langrefeq1077.gif)
added.
Deleting a row vector
![{z}](images/langref_langrefeq1077.gif)
from matrix
![{a}](images/langref_langrefeq35.gif)
corresponds to the rank-1 modification
![{c}^* = {c}- {z}^' {z}, { where } {c}^* = {a}^{*' {a}^*](images/langref_langrefeq1094.gif)
and the
![(m-1) x n](images/langref_langrefeq1095.gif)
matrix
![{a}^*](images/langref_langrefeq1096.gif)
contains
all rows of
![{a}](images/langref_langrefeq35.gif)
with the row
![{z}](images/langref_langrefeq1077.gif)
deleted.
Thus, you can also use the subroutines
RUPDT and
RDODT to update or downdate the Cholesky factor
![{{r}}](images/langref_langrefeq40.gif)
of a
positive definite crossproduct matrix
![{c}](images/langref_langrefeq1090.gif)
of
![{a}](images/langref_langrefeq35.gif)
.
The process of downdating an upper triangular matrix
![{{r}}](images/langref_langrefeq40.gif)
(and eventually a residual sum-of-squares
vector
ssq) is not always successful.
First of all, the downdated matrix
![{{r}}](images/langref_langrefeq40.gif)
could be rank deficient.
Even if the downdated matrix
![{{r}}](images/langref_langrefeq40.gif)
is of full rank,
the process of downdating can be ill-conditioned
and does not work well if the downdated matrix is
close (by rounding errors) to a rank-deficient one.
In these cases, the downdated matrix
![{{r}}](images/langref_langrefeq40.gif)
is not
unique and cannot be computed by subroutine RDODT.
If
![{{r}}](images/langref_langrefeq40.gif)
cannot be computed,
def returns 2, and the results
rup,
bup, and
sup return missing values.
The downdating of the residual sum-of-squares
vector
ssq can be a problem, too.
In practice, the downdate formula
![{ssq}_{new} = \sqrt{{{ssq}}_{old} - {{ssq}}_{dod}}](images/langref_langrefeq1097.gif)
cannot always be computed because, due to
rounding errors, the radicand can be negative.
In this case, the result vector sup
returns missing values, and def returns 1.
You can use various methods to compute the
columns
of the
matrix
that minimize the
linear
least squares problems with an
coefficient matrix
,
, rank
, and
right-hand-side vectors
(stored columnwise in the
matrix
).
The first of the following methods solves the
normal equations and cannot be applied to
the example with the
Hilbert matrix
since too much rounding error is introduced.
Therefore, use the following simple example:
a = { 1 3 ,
2 2 ,
3 1 };
b = { 1, 1, 1};
m = nrow(a);
n = ncol(a);
p = 1;
- Cholesky Decomposition of Crossproduct Matrix:
aa = a` * a; ab = a` * b;
r = root(aa);
x = trisolv(2,r,ab);
x = trisolv(1,r,x);
- QR Decomposition by Householder Transformations:
call qr(qtb,r,piv,lindep,a, ,b);
x = trisolv(1,r[,piv],qtb[1:n,]);
- Stepwise Update by Givens Rotations:
r = j(n,n,0.); qtb = j(n,p,0.); ssq = j(1,p,0.);
do i = 1 to m;
z = a[i,];
y = b[i,];
call rupdt(rup,bup,sup,r,z,qtb,y,ssq);
r = rup;
qtb = bup;
ssq = sup;
end;
x = trisolv(1,r,qtb);
Or equivalently:
r = j(n,n,0.);
qtb = j(n,p,0.);
ssq = j(1,p,0.);
call rupdt(rup,bup,sup,r,a,qtb,b,ssq);
x = trisolv(1,rup,bup);
- Singular Value Decomposition:
call svd(u,d,v,a);
d = diag(1 / d);
x = v * d * u` * b;
For the preceding
![3 x 2](images/langref_langrefeq1098.gif)
example matrix
![{a}](images/langref_langrefeq35.gif)
,
each method obtains the unique LS estimator:
ss = ssq(a * x - b);
print ss x;
To compute the (transposed) matrix
,
you can use the following specification:
r = shape(0,n,n);
y = i(m);
qt = shape(0,n,m);
call rupdt(rup,qtup,sup,r,a,qt,y);