Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
1/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
Organization (S):
EDF/IMA/MN
Manual of Reference
R6.03 booklet: Particular methods
Document: R6.03.01
Resolution of nonregular systems by
a method of decomposition in values
singular
Summary:
This document is devoted to the resolution of the systems of linear equations nonregular. Matrices taken
in account can be square noninvertible or rectangular.
After having recalled the theoretical framework of the solutions within the meaning of least squares, we concentrate the talk
on the method by decomposition in singular values which provides, on the one hand, a tool for diagnosis of the degree
of regularity of the system, and, in addition, a family of algorithms of resolution at the same time more general and more
stable that those drifting of the approach by the normal equations.
Lastly, we detail the algorithm implemented in Code_Aster which reabsorbs the equivalent functionality of
the bookstore Nag (F04JDF for version 12 and F04JDE for version 15) used for the modeling of
metallurgical behavior of steels [R4.04.01].
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
2/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
Contents
1 Introduction ............................................................................................................................................ 3
2 Solution of a rectangular linear system .......................................................................................... 4
2.1 Formalism of least squares .................................................................................................... 4
2.2 Existence of optima ..................................................................................................................... 5
2.3 Unicity of the optimum and row of the system ......................................................................................... 5
2.4 Solution within the meaning of least squares ............................................................................................ 6
3 singular Values ................................................................................................................................. 7
3.1 Decomposition in singular values ............................................................................................. 7
3.2 Row, image and core ..................................................................................................................... 8
3.3 Pseudo-opposite and solution within the meaning of least squares ............................................................... 9
4 Resolution of a rectangular linear system .................................................................................... 11
4.1 Method of the normal equations ................................................................................................. 11
4.2 Method by decomposition in singular values ...................................................................... 11
4.3 Comparison of the method of the normal equations to the method of decomposition in values
singular ..................................................................................................................................... 12
4.3.1 Conditioning .................................................................................................................. 12
4.3.2 Loss of precision ................................................................................................................ 12
4.3.3 Structure digs ................................................................................................................... 13
4.3.4 ............................................................................................................................ Conclusion 13
5 Algorithm SVD for the resolution of a linear équi or under-constrained system ................................ 13
5.1 Reduction of the problem and principle of the algorithm ....................................................................... 13
5.1.1 Reduction with the triangular form higher ....................................................................... 14
5.1.2 Reduction with the higher form bidiagonale ....................................................................... 14
5.1.3 Decomposition SVD of the bidiagonale higher ............................................................... 14
5.2 Reduction with the triangular form higher ................................................................................ 15
5.3 Reduction with the higher form bidiagonale ................................................................................ 17
5.4 Decomposition SVD of a higher bidiagonale ....................................................................... 19
5.4.1 Principle of the algorithm ........................................................................................................ 19
5.4.2 Implicit Diagonalisation of the normal matrix ................................................................... 20
5.4.3 Analyze decomposition ................................................................................................... 23
5.4.4 Organization of the algorithm ................................................................................................. 24
6 Bibliography ........................................................................................................................................ 24
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
3/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
1 Introduction
Being given a real matrix
With
of command
m N
×
and a vector
B
element of
¤
m
, we consider
the problem of the determination of a vector
X
element of R
N
who checks the following linear system:
Ax B
=
éq 1-1
It is well-known ([bib3] p. 9) that this system admits one, and only one solution, for all
B
element of
R
m
under the conditions necessary and sufficient that it is equi-constrained
(
)
m N
=
and that its matrix
With
that is to say regular. Also, the investigation of the under-constrained case
(
)
m N
and overstrained
(
)
m N
us
will confront with the one of the three following situations:
1) The linear system [éq 1-1] admits a solution and only one,
2) The linear system [éq 1-1] does not admit a solution,
3) The linear system [éq 1-1] admits an infinity of solutions.
In practice, situation 2) in general meets in the case of an overstrained system then
that the singular and under-constrained equi-constrained systems lead in general to situation 3).
To claim to solve a linear system of the type [éq 1-1], we should initially define what us
will call solution. This is the object of the paragraph 2 which is based mainly on the concept of
least squares and on differentiable optimization to define, some is the type of system,
a solution which is always single.
Paragraph 3 is devoted to the decomposition in singular values of the matrices (in summary
SVD: Been worth Singular Decomposition), which, not only constitutes a tool to diagnose which
of the three preceding situations corresponds to the studied linear system, but also provides a method
of determination of the solution defined in paragraph 2.
The method using decomposition SVD is presented at paragraph 4 and y is compared with
method of the normal equations.
Paragraph 5 details on the algebraic level the application of method SVD to the resolution of one
linear équi or under-constrained system such as it is implemented in Code_Aster.
In the following paragraphs, we will use the notations below:
·
X
and
()
X, y
for, respectively, the euclidian norm of the vector
X
and the scalar product
associated vectors
X
and
y
elements of
¤
m
or of
¤
N
,
·
M
T
for transposed of the matrix
M
,
·
Ker M
and
Imam
for, respectively, the core and the image of (the associated linear application
with) the matrix
M
,
·
X
for the orthogonal one of under space
X
of
R
m
or of
R
N
.
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
4/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
2
Solution of a rectangular linear system
In this paragraph we will define a concept of solution for the linear system [éq 1-1] which enjoys
properties of existence and unicity. The step proceeds in two times:
·
Initially, by an approach of the type least squares we build a problem
of differentiable and convex optimization (section 2.1) which admits always at least a solution
(section 2.2). Situation 2) of paragraph 1 is then eliminated,
·
Then, analyzing the property of unicity (section 2.3) to note that it is not always
guarantee we will impose an additional stress (section 2.4) on the solution
characterized in section 2.1 in order to restore unicity.
2.1
Formalism of least squares
The single solution of a linear system
Ax B
=
of square and regular matrix the minimum realizes of
quantity
Ay B
-
when
y
described
¤
N
. This property opens the channel to us which leads to a concept of
solution for a general linear system of the type [éq 1-1] which confers the same properties to him as
those of the particular case of the regular system. We will thus say point
X
of
¤
N
that it is solution of
system [éq 1-1] if it is solution of the problem of optimization:
Ax B
Ay B
y
- =
-
Min
U
N
éq 2.1-1
This approach is natural because it defines a solution of which the residue
R Ax B
=
-
is null in the case
where the second member is element of
Im A
and is of minimal standard in the contrary case, which
constitute best than one can wait.
To analyze the problem [éq 2.1-1], it is convenient to substitute the problem of optimization to him without
stresses are equivalent according to:
()
()
to find such
that
X
X
y
y
=
U
U
N
N
J
MinJ
éq 2.1-2
where
()
J.
is the functional calculus defined by:
()
J:
J
y
y
Ax B
=
-
U
N
1
2
The interest of the problem [éq 2.1-2] is due to the fact that the functional calculus
()
J.
check the following properties:
·
()
J.
is twice continuement differentiable:
()
()
(
)
DJ
:
DJ
X H
X H
In Ax A B, H
=
-
U
U
N
T
T
éq 2.1-3
()
()
()
()
(
)
D J
:
DJ
,
2
X
H, K
X H, K
In Ah K
×
=
U
U
U
N
N
T
éq 2.1-4
·
()
J.
is quadratic and convex.
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
5/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
Thus, the problem [éq 2.1-2] lies within the scope of the differentiable and convex optimization of kind
that we have the following results ([bib1] p. 156 and 146):
1) Convexity: any local optimum is in fact a total optimum, i.e. a solution of
[éq 2.1-2],
2) Differentiability: any local optimum checks the equation of Euler
()
DJ X
=
0
on
¤
N
who,
taking into account [éq 2.1-3], led to the characterization by the equations known as normal:
In Ax
With B
T
T
=
éq 2.1-5
2.2 Existence
optima
In [bib1] p. 171 one finds a demonstration of the existence of at least a solution to the equations
normals [éq 2.1-5]. This demonstration is based on arguments intended for the taking into account of
case of infinite dimension (theorem of projection on convex closed of a space of Hilbert).
Our case being definitely simpler, we give a demonstration of the result which only uses
simple algebraic arguments which, moreover, we will be useful in paragraph 3. To show that,
for all
B
element of
¤
m
, the normal equations [éq 2.1-5] admit a solution is equivalent to
establishment of inclusion
Im
Im
With
WITH A
T
T
. However, for any real matrix
M
of command
m N
×
us
let us have
(
)
Im
Ker
M
M
T
=
([bib3] p. 28). Also, the inclusion to be established which is equivalent to
(
)
(
)
Ker
Ker
With
WITH A
T
who is it even equivalent to
Ker
Ker
WITH A
With
T
. That is to say thus
X
WITH A
Ker
T
; then
Ax
With
Ker
T
, i.e.
(
)
Ax
With
Im
. Like
Ax
is also element of
Im A
, it can be only null what means that
X
With
Ker
and the demonstration completes.
This stage of the matter, we can say that any system of the type [éq 1-1] admits at least one
solution within the meaning of [éq 2.1-3] and all these solutions are characterized as solution within the meaning of
Cramer of the normal equations of [éq 2.1-5]. Situation 2) of paragraph 1 is eliminated.
Remain to eliminate situation 3), i.e. to guarantee unicity.
2.3
Unicity of the optimum and row of the system
It is clear that normal equations [éq 2.1-5], characterizing the optima which we seek,
admit a single solution under the condition necessary and sufficient that
WITH A
T
that is to say regular.
Like
WITH A
T
is always semi-definite positive, its inversibility is equivalent to its definite positivity, of
left that, taking into account the expression [éq 2.1-4] of the derived second of the functional calculus
()
J.
, us
let us find the well-known theorem of unicity of the optimum of the problem [éq 2.1-2] for a functional calculus
convex twice continuement differentiable ([bib1] HT 7.4-3 and 7.4-4).
In any general information, nothing prevents the matrix
WITH A
T
to be singular, the solution of the system [éq 1-1]
within the meaning of [éq 2.1-2] is thus not always single. We have nevertheless a criterion for
to detect this situation. With section 2.2 we established that
Im
Im
With
WITH A
T
T
and like
reciprocal inclusion is trivialement true, we can show the identity
Im
Im
With
WITH A
T
T
=
.
The introduction of the row
()
rg A
matrix
With
, the dimension of its space image, allows us then
of saying that a condition necessary and sufficient so that
WITH A
T
that is to say invertible is that
()
rg A A
T
N
=
what is equivalent to
()
rg A
=
N
because
() ()
()
rg
rg
rg
WITH A
With
With
T
T
=
=
.
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
6/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
The interest of this criterion is due to the fact that it limits the analysis to the only matrix
With
without it being necessary of
to form explicitly
WITH A
T
. This criterion also shows us that the normal equations associated one
strictly under-constrained linear system always admit an infinity of solutions. Indeed, it
row of a matrix is as equal to the number of independent columns as it has; also, for
that this row reaches the value
N
it is necessary that the columns of the matrix is of command at least
N
.
2.4
Solution within the meaning of least squares
We have just noted that the whole of the points which minimize the residue of the system [éq 1-1] is not
not necessarily reduced to only one point. To restore unicity we refine the concept of solution of
system [éq 1-1] of section 2.1 by defining the solution within the meaning of least squares like
the element of minimal standard of the whole of the points which minimize the residue. This solution
X
is
then characterized by:
{
}
X S
y
In Ax A B
X
y
T
T
y S
=
=
=
def
N
R;
Inf
and
This characterization is not satisfactory on the practical level because it asks for the resolution of one
problem of optimization under stresses. We will substitute another characterization more to him
adapted to the direction where it will lead (see section 4.2) to a procedure of calculation definitely simpler.
The unit
S
above is relocated of core of
WITH A
T
by any of the vectors solutions
normal equations [éq 2.1-5]. Also, the additional condition of minimization of the standard
be interpreted like a simple projection: the solution within the meaning of least squares of the system [éq 1-1]
is anything else only the projection of the origin of
¤
N
on the whole of the solutions of the equations
normals. Also, we can characterize it like the point of intersection between the unit
S
and
the orthogonal one of the core of
WITH A
T
.
The definition of a solution to the system [éq 1-1] can then be summarized as follows:
(
)
X
Ax B
In Ax
With B
X
WITH A
is solution of
=
=
T
T
T
Ker
éq 2.4-1
The first condition makes
X
a vector of minimal residue while the second selects,
among the vectors of minimal residue, that of minimal standard.
The definition [éq 2.4-1] is a conventional generalization of the concept of solution of a system
equi-constrained regular and confers on any system of the type [éq 1-1] a solution and only one.
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
7/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
3 Values
singular
In this paragraph we have some useful results for the design of a method of
operational resolution of the system [éq 1-1]. These results derive from the concept of singular values
(section 3.1) and allow to build a base of the core and a base of the image of the matrix of
system (section 3.2) from which it is possible to give a direction, adapted to the calculation of
solution within the meaning of least squares, contrary to an unspecified matrix (section 3.3).
3.1
Decomposition in singular values
Let us start by pointing out the definition of the singular values. One calls singular values of one
real matrix
With
of command
m N
×
square roots of the eigenvalues of the square matrix
WITH A
T
of command
N
who, let us recall it, is semi-definite positive.
The concept of diagonalisation of the square matrices (when they are diagonalisables) spreads with
rectangular matrices (without restriction) by the concept of decomposition (or factorization) in values
singular.
For all real matrices
With
of command
m N
×
, there are two unit square matrices
Q
and
P
of a respective nature
m
and
N
such as:
WITH Q P
=
T
éq 3.1-1
where
is a matrix of command
m N
×
the structure is schematized below:
µ
µ
µ
1
2
!
N
0
=
if
m N
µ
µ
µ
1
2
!
N
0
=
if
m N
>
µ
I
are the singular values of
With
that we suppose ordered by descending order:
µ
µ
µ
1
2
“
N
One can find a demonstration of this result in [bib1] p.10 for the équi-constrained case and in
[bib3] p.73 for the overstrained case, the under-constrained case results some then by transposition.
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
8/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
Factorization SVD [éq 3.1-1] of
With
give
WITH A P
P
T
T
T
=
and
AA
Q
Q
T
T
T
=
so that,
T
and
T
being diagonal square matrices, the matrix
P
is consisted of the vectors
clean orthonormalized of the matrix
WITH A
T
while the matrix
Q
is consisted of the vectors
clean orthonormalized of the matrix
AA
T
.
3.2
Row, image and core
The paragraph [§2] showed the fundamental part which the row of the matrix plays
With
and the core of
stamp
WITH A
T
for the resolution of a nonregular linear system of the type [éq 1-1]. We will see
now how factorization [éq 3.1-1] can be used to determine this row like one
base
Ker A A
T
.
That is to say
R
the index of the smallest nonnull singular value. Factorization [éq 3.1-1] is written too
Q AP
T
=
where the taking into account of the null singular values makes it possible to specify the decomposition
in block of
:
R
0
=
if
m N
0
0
R
0
=
if
m N
>
0
0
0
where
(
)
R
R
=
Diag
,
,
µ µ
µ
1
2 “
is the diagonal matrix of command
R
nonnull singular values
in the ascending order.
Since matrices
Q
and
P
are regular, the matrices
With
and
are equivalent so that them
respective core and image coincide. We thus deduce from it that:
·
The row of
With
coincide with the number of nonnull singular values:
rg A
=
R
·
The vectors columns of
P
of index
R
+
1
with
N
form a base of
Ker A
·
The vectors columns of
Q
corresponding to the nonnull singular values form one
base
Im A
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
9/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
In addition, with section 2.3 we saw that
Im
Im
With
WITH A
T
T
=
. Identity
(
)
Im
Ker
M
M
T
=
us gives then
Ker
Ker
With
WITH A
=
T
so that the second condition of the definition [éq 2.4-1] is
simply realized by any vector which is expressed like a linear combination of the vectors
columns of
P
corresponding to the nonnull singular values.
3.3
Pseudo-opposite and solution within the meaning of least squares
Another application of the decomposition in singular values consists of the concept of
pseudo-opposite (or opposite Moore-Penrose) which generalizes the usual concept of reverse of a matrix
square regular with the rectangular matrices on the one hand, and the singular square matrices of other
leaves.
First of all, theopposite one of a matrix
decomposition in singular values [éq 3.1-1]
is defined by:
-
R 1
0
=
+
if
m N
0
0
-
R 1
0
=
+
if
m N
>
0
0
0
where
R
1
2
R
-
=
1
1
1
1
Diag
,
,
µ µ
µ
#
is the reverse with the usual direction of
R
.
This being, we use the decomposition [éq 3.1-1] matrix
With
to define its pseudo-opposite
With
+
by:
With
P
Q
+
+
=
T
éq 3.3-1
In the same way, of the decomposition [éq 3.1-1] of the matrix
With
we fire
WITH A P
P
T
T
T
=
, of kind
that theopposite one
()
WITH A
T
+
matrix
WITH A
T
is defined by:
()
()
WITH A
P
P
T
T
T
+
+
+
=
éq 3.3-2
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
10/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
We are able now to provide a simple interpretation of the solution within the meaning of
least squares defined by [éq 2.4-1].
The restriction on
(
)
Ker A A
T
linear application associated the matrix
WITH A
T
one defines
isomorphism of
(
)
Ker A A
T
on
Im A A
T
. Like, on the one hand
(
)
(
)
Ker
Ker
Im
WITH A
With
With
T
T
=
=
, and, in addition,
Im
Im
WITH A
With
T
T
=
, this restriction is in fact
an automorphism of
(
)
Ker A A
T
. In the base of
(
)
Ker A A
T
constituted by
R
first
columns of the matrix
P
, this automorphism is represented by the matrix
R
2
. Also, sound
reciprocal automorphism is represented there by the matrix
R
-
2
. Extension to
¤
N
of this
automorphism is then represented, in the base associated with the matrix
P
, by the matrix
()
()
T
T
T
+
+
=
, and thus, in the canonical base, by the matrix
()
WITH A
T
+
.
It follows that:
·
We find the fact that, for all
B
element of
¤
m
, there is a single vector
(
)
X
WITH A
Ker
T
solution of
In Ax A B
T
T
=
, that is to say the existence and the unicity of the solution with
feel least squares [éq 2.4-1) of the system
Ax B
=
,
·
This single solution is given by:
()
X
WITH A A B
=
+
T
T
éq 3.3-3
Theopposite one of a matrix is defined starting from decomposition SVD of this matrix. Like
decomposition SVD is not single, the pseudo-opposite matrix is not single. On the other hand, of
point of view of the linear applications associated the matrices, the pseudo-opposite application is single.
All the matrices pseudo-opposite associated with various decompositions SVD with a matrix
data are not whereas representing matric particular which expresses this application
pseudo-opposite relative at the bases induced by the orthogonal matrices of the decompositions
SVD. Also, the expression [éq 3.3-3] has a direction: it defines a vector of which
X
represent them
components compared to the base of arrival (matrix
P
) of decomposition SVD.
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
11/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
4
Resolution of a rectangular linear system
The two methods of resolution of the system [éq 1-1] which we present at sections 4.1 (method
normal equations) and 4.2 (decomposition in singular values) aim to the resolution of
normal equations [éq 2.1-5]. These two methods are characterized not only by the choice of
algorithms which they implement (inversion against pseudo-inversion), but also by their degree
of general information and by their numerical properties which are compared with section 4.3
4.1
Method of the normal equations
The resolution of the system
Ax B
=
by the method of the normal equations consists in calculating
solution within the meaning of least squares [éq 2.4-1] in a “direct” way, i.e. while using directly
the relation
()
X
WITH A
With B
=
-
T
T
1
. For that, it is initially a question of calculating
WITH A
T
and
With B
T
, then to solve
the system obtained either by iterative method or by a factorization of
WITH A
T
.
We can notice right now that this method is limited to the matrices
WITH A
T
regular,
what limits its applicability to the system [éq 1-1] whose matrix is of full row (see
section 2.3). In particular, the method of the normal equations cannot treat nor the systems
strictly under-constrained, nor singular équi-constrained systems (see section 2.4).
4.2
Method by decomposition in singular values
We saw with section 3.3 that the solution of the system
Ax B
=
within the meaning of least squares
defined by [éq 2.4-1] can be characterized by the relation
()
X
WITH A A B
=
+
T
T
[éq 3.3-3]. Method
of resolution of the system based on this property is known as method by decomposition in values
singular because it builds theopposite one [éq 3.3-2] of
WITH A
T
via decomposition SVD
[éq 3.1-1] of the matrix
With
.
As any matrix can be broken up into singular values, it follows that any system of the type
[éq 1-1] can be solved within the meaning of [éq 2.4-1] by this method which thus presents, at least,
the advantage of the general information compared to the method of the normal equations.
It is not all. Method by decomposition in singular values, contrary to the method
normal equations, does not require the explicit construction of the matrix
WITH A
T
and of the vector
With B
T
(we will see with section 4.3 the interest on the numerical level of this property). Indeed, it is
easy to check that the matrix
singular values of factorization [éq 3.1-1] satisfied with
identity
()
+
+
+
+
T
T
+
+
=
, so that, prémultipliant
With
T
by theopposite one of
WITH A
T
and
taking account of factorization [éq 3.1-1], we obtain
()
()
WITH A A
P
P P
Q
T
T
T
T T
T
T
+
+
+
=
, which, by orthogonality of
P
us gives
()
WITH A A
With
T
T
+
+
=
. Consequently, the characterizations [éq 3.3-3] and [éq 2.4-1] of the sought solution are
equivalent to the characterization:
X
Ax B
X With B
is solution of
= =
+
éq 4.2-1
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
12/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
4.3 Comparison of the method of the normal equations to the method
of decomposition in singular values
In the two preceding sections, we come to note, that algebraically, the method of
decomposition in singular values is more general and simpler than the method of the equations
normals. We now will note, while following [bib3] p. 336, that it is also higher to him
on the numerical level. This superiority is expressed on the one hand, in term of stability not only of
resolution, but also of construction of the problem and, in addition, on a less critical level, in
term of adaptation to the processing of the hollow matrices.
4.3.1 Conditioning
The conditioning of a matrix
With
of command
m N
×
is defined like the report/ratio of its values
singular extremes and nonnull:
()
cond A
= µµ
1
R
where
R
is the row of the matrix
With
.
Results presented in [bib3] p.184, using the normal equations as a tool for analysis and
not like a computational tool, show that the disturbance of the solution of the problem of optimization
[éq 2.1-1] due to the round-off errors can be proportional to
()
cond A
2
. But results
conventional of the analysis of stability of the solution of a linear system compared to these same
errors show a proportionality with the number of conditioning of the matrix. So that in
case of a direct resolution of the normal equations, we obtain an error always proportional
with
()
()
cond
cond
WITH A
With
T
=
2
, which is worse than
()
cond A
.
The method of resolution by decomposition in singular values uses only transformations
orthogonal (see paragraph 4), so that it does not modify the initial conditioning of
problem [éq 1-1] and is thus, from this point of view, more attractive than the method of the equations
normals.
4.3.2 Loss of precision
We have just seen that the round-off errors lead to a degradation of the solution more
sensitive when it is calculated via the normal equations rather than by a decomposition in
singular values. The following example, drawn from [bib 2], shows that construction even system
[éq 2.1-5] of the normal equations is disturbed by the round-off errors.
That is to say thus the following matrix:
With
=
1 1
0
0
leading to
WITH A
T
= +
+
1
1
1
1
2
2
whose singular values are
µ
1
2
=
+
2
and
µ
2
=
, so that the row of
With
is 2 as soon as
0
. If
check
2
Mach
<
<
where
Mach
is the precision machine, then all the coefficients of
WITH A
T
will be calculated with value 1 and the calculated singular values will be, at best,
µ
1
=
2
and
µ
2
=
0
. It follows that the numerical row, calculated by the normal equations, will be 1, whereas that
calculated by a decomposition SVD of the matrix
With
would be equal to 2
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
13/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
4.3.3 Structure
hollow
With a less level, the construction of the matrix of the normal equations induced a filling of
associated system that the method using the decomposition in singular values avoids.
4.3.4 Conclusion
The following table summarizes the discussion of the preceding sub-sections:
General information
Conditioning
Loss of precision with
construction of the problem
filling
Equations
normals
systems of full row
()
cond A
2
possible
yes
SVD
any system
()
cond A
impossible
not
5
Algorithm SVD for the resolution of a linear system
équi or under-constrained
In this paragraph, we detail the method of resolution of the not-regular systems put in
work in Code_Aster. This method applies to the under-constrained or équi-constrained system
singular and provides the solution within the meaning of least squares [éq 2.4-1].
The calculation of a decomposition SVD of
With
is equivalent to the calculation of the spectrum of the normal matrix
associated
WITH A
T
. Also, it can be obtained only with the convergence of an iterative process.
Section 5.1 exposes the principle of the algorithm and shows in particular how the application of
two orthogonal transformations makes it possible to reduce the problem to the simple search of
decomposition SVD of a higher matrix bidiagonale. Sections 5.2 and 5.3 are devoted to
the algorithmic one of these reductions. Section 5.4 presents the algorithm of decomposition SVD
matrix bidiagonale.
The algorithms will be described with the convention of notation in which:
·
()
R,
I J
indicate the rotation of Givens of the plan
()
I J
,
and of angle
,
·
()
With
K
indicate reiterated index
K
of a matric iteration and
()
With
K L,
reiterated
L
of an iteration
intern with reiterated
()
With
K
.
5.1
Reduction of the problem and principle of the algorithm
In this section, we present the algorithm of resolution of a linear sytème équi or
under-constrained by method SVD.
We reduce the problem to the search of decomposition SVD of a matrix bidiagonale like
in [bib2] but we carry out the reduction in another way that that proposed in [bib2]: us
let us start by reducing the matrix to a higher triangular form, then, we reduce this
triangular with a higher form bidiagonale. These two reductions are carried out by
orthogonal transformations.
The arithmetic operations of decomposition SVD are connected as follows:
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
14/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
5.1.1 Reduction with the higher triangular form
[
]
With
U
P
With
U P
1
1
=
<
=
=
0
T
T
m N
m N
if
if
éq 5.1-1
where
P
1
is an orthogonal matrix of command
N
and
U
a higher triangular matrix of command
m
.
5.1.2 Reduction with the higher form bidiagonale
U
Q LP
2
2
=
T
éq 5.1-2
where
Q
2
and
P
2
are two orthogonal matrices of command
m
and
B
a higher matrix bidiagonale
of command
m
.
5.1.3 Decomposition SVD of the higher bidiagonale
B Q
P
3
3
=
T
éq 5.1-3
where
Q
3
and
P
3
are two orthogonal matrices of command
m
and
a diagonal matrix of command
m
form:
µ
µ
1
2
!
0
0
!
0
0
=
R
0
0
=
0
Combining the relations [éq 5.1-1], [éq 5.1-2] and [éq 5.1-3], we obtain a decomposition SVD of
the matrix
With
:
WITH Q Q
=
2
3
R
0 0
0
0
P
1T
P P
3
2
T
T
0
0
I
éq 5.1-4
The solution within the meaning of least squares [éq 2.4-1] of the system [éq 1-1] is then obtained by
the application of pseudo-opposite [éq 3.3-1] of
With
deduced from the decomposition in singular values
[éq 5.1-4]. We thus obtain:
X P P P
Q Q B
2 3
3
2
=
+
1
0
T
T
éq 5.1-5
The algorithm proposed thus consists of the sequence of factorizations [éq 5.1-1], [éq 5.1-2] and
[éq 5.1-3] before the application of the relation [éq 5.1-5].
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
15/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
5.2
Reduction with the higher triangular form
Starting from a matrix
With
of command
m N
×
for
m N
a higher triangular matrix is determined
U
of command
m
and an orthogonal matrix
P
of command
N
such as:
[
]
With
U
P
WITH U P
=
<
=
=
0
T
T
m N
m N
if
if
Algorithmiquement, factorization uses a method of elimination which is interpreted, like
construction of a succession of matrices
()
With
K
by:
()
(
)
() ()
With
With
With
WITH P
m
K
K
K
K
m m
=
=
=
-
-
1
1
1
for
,
,
#
where each current matrix
()
With
K
have the structure schematized below:
Coefficients
cancelled in
phase 2
(line K)
Pivot (coefficient (K, K))
Column pivot
(column K)
Coefficients
not changing
more
coefficients
null
With
(K)
=
Coefficients
cancelled in
phase 1
(line K)
Coefficients
()
has
I J
K
,
matrices
()
With
K
iteration thus check:
()
has
K
I m
m
J N
K
I m
J K
K
J I m
I J
K
,
=
+
+
+
+
0
1
1
1
1
1
if
and
and
éq 5.2-1
so that at the end of the recursion, we will have:
()
()
U WITH
P
P
=
=
=
-
1
1
1
and
K
K m m
,
,
#
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
16/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
The problem is thus reduced to the safeguarding of the structure [éq 5.2-1] at the time of the passage of
()
With
K
with
(
)
With
K
-
1
by a transformation
()
P
K
who must be orthogonal. The problem of orthogonality is regulated
by choosing the transformation like a product of rotations of Givens and the problem of
safeguarding of the structure is solved by carrying out this product in an order which does not destroy them
zeros created.
Taking account of the rectangular structure of the matrix
()
With
K
, we build reiterated
(
)
With
K
-
1
in
two phases:
·
Phase 1 cancels successively the coefficients
(
)
has
K J
K,
-
1
corresponding to the columns
J K
K
= -
-
1
2
1
,
,
#
, which results in:
(
)
()
(
)
(
)
()
(
)
With
With
With
With
K
K
K
K
J
K
J
J K
T
K J
K
K
J
- -
- -
-
=
=
= -
-
1
1
1
1
1
1
2
1
,
,
,
R,
,
,
for
#
·
Phase 2 cancels successively the coefficients
(
)
has
K J
K,
-
1
corresponding to the columns
J N N
m
=
-
+
,
,
,
1
1
#
, which results in the recursion:
(
)
(
)
(
)
(
)
()
(
)
With
With
With
With
K
N
K
K
J
K
J
J K
T
K J
N N
K
J
-
-
- -
-
=
=
=
-
+
1
1 0
1
1
1
1
1
,
,
,
,
R,
,
,
for
#
The angle
()
J K
rotation of Givens of the plan
()
K J
,
is selected to cancel the coefficient in position
()
K J
,
of
(
)
has
K
J
-
1,
. The application of each rotation thus modifies only the columns
K
and
J
what
does not destroy the null coefficients produced by the preceding stages. We note that
column
K
play a particular part (that of pivot) because it only is systematically modified by
each rotation whereas the other columns are modified only by the rotation which cancels them
coefficient with the line
K
.
With the exit of these recursions, we have
(
)
(
)
With
With
K
K
K
-
-
=
1
1,
. The matrix
()
P
K
is then given by:
()
()
(
)
()
(
)
P
K
J m
J N
J K
J
J K
J K
K J
K J
=
= +
=
=
= -
1
1
1
R,
R,
so that the matrix
P
is worth:
()
(
)
()
(
)
P
=
=
=
= +
=
=
= -
K
K m
J m
J N
J K
J
J K
J K
K J
K J
1
1
1
1
R,
R,
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
17/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
5.3
Reduction with the higher form bidiagonale
To reduce a higher triangular square matrix
With
of command
m
with the higher form bidiagonale
consist in finding two matrices orthogonal
P
and
Q
and a higher matrix bidiagonale
B
,
all three of command
m
, such as:
WITH QBP
=
T
Algorithmiquement, factorization proceeds like that of the preceding section by using one
method of elimination which is interpreted algebraically like the construction of a succession of matrices
()
With
K
by:
()
(
)
()
() ()
With
With
With
Q
WITH P
1
1
1 2
2
=
=
=
-
+
K
K
K
K
T
K
m
for
,
#
where each current matrix
()
With
K
have the diagonal structure per block following:
·
The diagonal block higher (indices of line and column varying of 1 than
K
-
1
) is one
stamp bidiagonale higher command
K
-
1
,
·
The lower diagonal block (indices of line and column varying from
K
with
m
) is a matrix
triangular higher of command
m K
-
.
Coefficients
()
has
I J
K
,
matrices
()
With
K
iteration thus check:
()
has
I K
J m
I K
J
K I m
J I
I J
K
I
I
,
=
-
+
-
<
<
0
1
1
2
1
1
1
if
and
and
and
éq 5.3-1
so that at the end of the recursion, we will have:
(
)
()
()
B WITH
Q
Q
P
P
=
=
=
+
=
=
=
=
m
K
K m
K
K
K m
K
1
1
1
,
and
As for the factorization of the preceding section, the problem is reduced to safeguarding
structure [éq 5.3-1] at the time of the passage of
()
With
K
with
(
)
With
K
+
1
. The orthogonality of the transformations
()
Q
K
and
()
P
K
is obtained by building them like product of rotations of Givens and the problem of
safeguarding of the structure is solved by carrying out these products in an order which does not destroy them
zeros created by the preceding stages.
The algorithm thus cancels successively the coefficient
(
)
K J
,
+
1
for
J m
m
K
= -
-
+
1
2
2
,
,
#
by
the application on the right of a rotation of Givens of the plan
(
)
J J
,
+
1
. This rotation modifies only them
columns
J
and
J
+
1
, which creates a parasitic coefficient in position
(
)
J
J
+
1,
. This parasitic coefficient
is then eliminated by the application on the left from transposed of a rotation of Givens in the plan
(
)
J J
,
+
1
.
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
18/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
The process of passage of
()
With
K
with
(
)
With
K
+
1
can then be formalized by:
(
)
()
(
)
(
)
()
(
)
(
)
()
(
)
(
)
(
)
(
)
With
With
With
With
With
With
With
With
K
m
K
K
J
K
J
J K
K
J
J K
T
K
J
K
K
K
J J
J J
J m
m
K
+
-
+ =
+
+ =
=
+ =
+
+
=
=
+
=
+
= -
-
+
=
1
1
1
1 2
1
1
1
1 2
1
1 2
1
1
1
1
1
2
1
,
,
/
,
,
/
,
/
,
R,
,
R,
,
,
,
for
#
where angles
()
J K
and
()
J K
=
1 2
/
are selected to cancel the coefficient in position respectively
(
)
K J
,
+
1
of
(
)
With
K
J
+
1,
and the coefficient in position
(
)
J
J
+
1,
of
(
)
With
K
J
+ =
1
1 2
,
/
.
The structure of the matrices
(
)
With
K
J
+
1,
is illustrated in the following figure:
Line
K
Bi-diagonal
Coefficient to be cancelled by rotation
P
current (column
j+1
)
Coefficient cancelled by
rotations
P
following
Coefficient cancelled by
rotations
P
the preceding ones
With
(,)
K J
=
Coefficient created by
rotation
P
current and
eliminated by rotation
Q
current
With the exit of this recursion, matrices
()
P
K
and
()
Q
K
are given by:
()
()
(
)
()
()
(
)
P
Q
K
J m
J K
J K
K
J m
J K
J K
J J
K J
=
+
=
= -
= +
= -
= +
=
1
1
1
1
1 2
1
R,
,
R,
/
and
so that matrices
P
and
Q
are worth:
()
(
)
()
(
)
P
Q
=
+
=
+
=
= -
= -
= +
=
= -
= -
= +
=
K
K m
J m
J K
J K
K
K m
J m
J K
J K
J J
J J
1
2
1
1
1
2
1
1
1 2
1
1
R,
,
R,
,
/
and
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
19/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
5.4
Decomposition SVD of a higher bidiagonale
We present an algorithm of construction of decomposition SVD of a matrix bidiagonale
higher
With
of command
m
. The algorithm thus builds two orthogonal matrices
Q
and
P
and one
stamp diagonal
D
such as:
WITH QDP
=
T
The algorithm is drawn from [bib2].
5.4.1 Principle of the algorithm
The calculation algorithm of decomposition SVD diagonalise repeatedly the matrix
With
by means of
recursion:
()
(
)
()
() ()
With
With
With
Q
WITH P
1
1
1 2
=
=
=
+
K
K
K
K
T
K
for
, #
éq 5.4.1-1
where matrices
()
Q
K
and
()
P
K
are orthogonal and the matrices
()
With
K
are bidiagonales higher.
With convergence, we will have:
()
()
()
D WITH
P
P
Q
Q
=
=
=
=
=
=
=
,
K
K
K
K
K
K
1
1
and
The idea of the iteration consists with:
·
To choose
()
P
K
to make converge the algorithm
QR
applied to the diagonalisation of the matrix
(known as normal)
WITH A
T
without forming it explicitly. Indeed, the matrix
P
decomposition
SVD of
With
is anything else only the matrix of the clean vectors of
WITH A
T
,
·
To choose
()
Q
K
to preserve the higher structure bidiagonale the reiterated successive ones.
As in the case of the factorizations presented at sections 5.2 and 5.3, the matrices
()
Q
K
and
()
P
K
are built like product of rotations of Givens. The passage of
()
With
K
with
(
)
With
K
+
1
is then
realized by:
(
)
() ()
()
[
]
()
() ()
()
[
]
With
Q
Q
Q
With
P
P
P
K
K
K
K m T
K
K
K
K m
+
=
=
1
2
3
2
3
,
,
,
,
,
,
#
#
éq 5.4.1-2
where them
()
Q
K I,
and
()
P
K I,
are two rotations of the plan
(
)
I
I
-
1,
of respective angle
J
and
I
:
()
()
(
)
()
()
(
)
Q
P
K I
I K
K I
I K
I
I
I
I
,
,
R
,
R
,
=
-
=
-
1
1
and
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
20/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
Rotations are alternatively applied on the right then on the left so that
(
)
With
K
+
1
preserve the higher structure bidiagonale of
()
With
K
. With this intention:
The angle
()
2k
is, for the moment, arbitrarily selected; the application of rotation
()
P
K, 2
create
then a coefficient in position
()
2 1
,
,
The angle
()
2 K
is selected so that the application of rotation
()
Q
K, 2
cancel the coefficient in position
()
2 1
,
, which creates a coefficient not no one in position
()
1 3
,
,
The angle
()
3k
is selected so that the application of rotation
()
P
K, 3
cancel the coefficient in position
()
1 3
,
, which creates a coefficient not no one in position
()
3 2
,
,
$
The angle
()
mk
-
1
is selected so that the application of rotation
(
)
Q
K m
,
-
1
cancel the coefficient in
position
(
)
m
m
-
-
1
2
,
, which creates a coefficient not no one in position
(
)
m
m
-
2,
,
The angle
()
mk
is selected so that the application of rotation
()
P
K m
,
cancel the coefficient in position
(
)
m
m
-
2,
, which creates a coefficient not no one in position
(
)
m m
,
-
1
,
The angle
()
mk
is selected so that the application of rotation
()
Q
K m
,
cancel the coefficient in
position
(
)
m m
,
-
1
, and stamps it
(
)
With
K
+
1
that is to say bidiagonale higher.
For any value of the angle, this process ensures maintains it structure bidiagonale higher than
reiterated [éq 5.4.1-2]. We will see now how it is possible to choose this angle to make
to converge the iteration [éq 5.4.1-1].
5.4.2 Implicit Diagonalisation of the normal matrix
The algorithm of the preceding sub-section leaves unspecified the angle
()
2k
the first rotation of
()
P
K
. We will raise this indetermination in order to make matrix
()
P
K
the matrix
orthogonal of a pitch
QR
, with spectral shift, applied to the diagonalisation of the normal matrix
MR. A A
T
=
.
With iteration SVD [éq 5.4.1-1] of the matrix
With
, we associate an iteration on the normal matrix
MR. A A
T
=
:
(
)
(
)
(
)
()
() ()
M
With
With
P
MR. P
K
K
K
K
K
K
T
T
+
+
+
=
=
1
1
1
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
21/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
Iteration
QR
for the diagonalisation of the normal matrix
The transformation
QR
, with spectral shift
K
, applied to
()
M
K
is written:
To factorize
()
M
I
K
K
-
in the form
()
M
I P R
K
K
-
=
To build
(
)
M
K
+
1
by
(
)
M
P R
I
K
K
+
=
+
1
where
P
and
R
are two matrices respectively orthogonal and triangular higher. Matrices
()
M
K
and
(
)
M
K
+
1
are thus tridiagonales and similar:
(
)
()
M
P MR. P
K
T
K
+
=
1
From the practical point of view, the matrix
P
is presented in the form of a product of rotations of Givens:
(
) (
) (
)
P
T
N
N
N
N
N
N
=
-
-
-
-
R
,
R
,
,
R,
1
2
1
1 2
1
2
#
Angles
K
are selected so that the application on the left of
(
)
R
,
K
K
K
-
1
with the matrix
(
)
()
(
)
1
1
2
2
1
1 11
= - -
-
-
K
K
K
K
,
,
R
,
#
M
I
cancel the coefficient of position
(
)
K K
,
-
1
in the matrix
result.
Francis showed that the passage of
()
M
K
with
(
)
M
K
+
1
do not require the formation clarifies
stamp
()
M
I
K
K
-
: the shift can be carried out implicitly. The theorem is stated like
follows:
Theorem (Francis): That is to say
X
an orthogonal matrix whose first column coincides with that
of
P
. Under the assumptions:
1)
(
)
()
M
X MR. X
K
T
K
+
=
1
2)
(
)
M
K
+
1
is tridiagonale,
3) Under-diagonal elements of
()
M
K
are all nonnull (irreducibility of
()
M
K
),
one has
(
)
(
)
M
DM
D
K
K
+
+
=
1
1
where
D
is a diagonal matrix of diagonal coefficients all equal to
±
1
.
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
22/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
Application to algorithm SVD
Consequently, the choice of the angle
()
2k
the first rotation of iteration SVD [éq 5.4.1-1] consists of
()
2
2
K
= -
. Thus
()
(
)
(
)
R,
R,
1 2
1 2
2
2
K
T
=
so that the first column of
()
P
K
coincide
with the first column of
P
. Therefore, if all under-diagonal elements of
()
M
K
are nonnull,
then, matrices
()
P
K
and
P
are identified (with a multiplicative factor
±
1
close to the columns) and
iteration SVD [éq 5.4.1-1] is equivalent to the application of the transformation
QR
, with shift, with
stamp
()
M
K
.
Choice of the spectral shift
The shift is usually selected like the eigenvalue of the lower miner of command two of
()
With
K
nearest to
()
has
m m
K,
. This choice ensures a total convergence which, generally, is
cubic.
Applicability of the theorem of Francis and phenomenon of decomposition
The use of the theorem of Francis supposes nonthe nullity of all the under-diagonal coefficients of
matrices
()
M
K
, which of anything is not guaranteed. Moreover, within the framework of the method
QR
of
diagonalisation of a symmetrical matrix tridiagonale, appearance of null under-diagonal coefficients
is:
·
Desirable: the null coefficients uncouple the diagonal blocks which they frame, which
bring back the diagonalisation of the complete matrix to the diagonalisation of its diagonal blocks
(this phenomenon is often called “decomposition”),
·
Inevitable: the convergence of the algorithm towards an eigenvalue is interpreted algebraically
like the appearance of a preceding diagonal block of command 1.
In the following sub-section, we will see which processing it is appropriate to adopt in the presence of one
decomposition.
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
23/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
5.4.3 Analyze decomposition
The analysis of decomposition relates to each reiterated taken independently of the others, also, us
will not use the superscript
()
K
.
That is to say
D D
D
m
1
2
,
,
#
and
E E
E
m
2
3
,
#
respectively diagonal and on-diagonal elements of
With
.
Under-diagonal elements of the normal matrix
MR. A A
=
T
are then given by:
m
D E
I
m
I
I
I L
+
+
=
=
-
1
1
1 2
1
,
,
for
#
Let us suppose, to simplify, that only the coefficient
m
L
L
-
1,
is null. The matrix
M
present one then
structure of two diagonal blocks whose meeting of the respective spectra gives the spectrum of
M
.
This decomposition takes place is for
E
1
0
=
maybe for
E
D
L
1
1
0
0
=
-
and
.
The case
E
1
0
=
do not raise any difficulty. The matrix
With
then have a diagonal structure of two
blocks which provide each one a part complementary to decomposition SVD of
With
. Each block
being bidiagonal higher, without null on-diagonal coefficients, its decomposition SVD is calculable
by the iteration [éq 5.4.1-1].
The case
E
D
L
1
1
0
0
=
-
and
is more delicate. Indeed, the iteration [éq 5.4.1-1] cannot be applied nor
with the matrix
With
, to violate the assumptions of the theorem of Francis, nor with none under matrices
of
With
, to ensure the structure bidiagonale reiterated. This problem is circumvented by
postmultiplication of
With
by a series of rotations of Givens in the successive plans
(
) (
) (
)
L
L L
L
L
m
-
-
+
-
1
1
1
1
,
,
,
,
#
:
·
The rotation of the plan
(
)
L
L
-
1,
cancel the coefficient
(
)
L
L
-
1,
and creates a coefficient in position
(
)
L
L
-
+
1
1
,
,
·
The rotation of the plan
(
)
L
L
-
+
1
1
,
cancel the coefficient
(
)
L
L
-
+
1
1
,
and creates a coefficient in
position
(
)
L
L
-
+
1
2
,
,
·
The rotation of the plan
(
)
L
L
-
+
1
1
,
cancel the coefficient
(
)
L
L
-
+
1
1
,
and creates a coefficient in
position
(
)
L
L
-
+
1
2
,
,
$
·
The rotation of the plan
(
)
L
m
-
1,
cancel the coefficient
(
)
L
m
-
1,
and does not create a coefficient.
So that the matrix produced by this process has the same structure as that corresponding
with the case
E
1
0
=
.
Code_Aster
®
Version
3.0
Titrate:
Resolution of nonregular systems by a method of decomposition
Date:
26/06/95
Author (S):
B. QUINNEZ, R. MICHEL
Key:
R6.03.01-A
Page:
24/24
Manual of Reference
R6.03 booklet: Particular methods
HI-75/95/032/A
5.4.4 Organization of the algorithm
The algorithm isolates successively each singular value, also, it exists an index
K
p
such as reiterated
()
With
K
p
breaks up into two diagonal blocks:
()
With
K
p
=
()
B
K
p
0
0
()
D
K
p
where
()
B
K
p
is a higher matrix bidiagonale of command
p
and
()
D
K
p
is a diagonal matrix
of command
m p
- +
1
gathering on its diagonal the found singular values.
From this reiterated, the algorithm applies the iteration of sub-section 5.4.1 to the submatrix
()
B
K
p
until the cancellation of the coefficient in position
(
)
p
p
-
1,
, signal of the convergence of
p
ième
value
singular. Each pitch of the internal iteration, thus defined, is organized as follows:
·
Analyze decomposition,
·
If the found decomposition corresponds to a diagonal element no one, a series of
additional rotations is applied to find the structure of generated decomposition
by a on-diagonal element no one. These rotations are built according to the method presented
with sub-section 5.4.3,
·
If the coefficient in position
(
)
p
p
-
1,
then does not produce decomposition the submatrix
()
B
K
p
is the object of a pitch of the iteration [éq 5.4.1-1] where, in accordance with the analysis of under
section 5.4.2, the angle of the first rotation is selected so that this pitch is equivalent to
the application of a transformation
QR
, with implicit spectral shift, on the matrix
associated normal.
The complete convergence of the iteration is then obtained with the index
K
m
for which the submatrix
()
D
K
m
is of command
m
.
Of course, in practice, a coefficient is regarded as null as soon as it is lower, in value
absolute, with a certain tolerance. The tolerance generally used for the problems of values
singular is selected like the product of the precision machine by
With
1
. Let us notice that in
case of a decomposition produced by a diagonal element no one with the tolerance chosen, the application of
series of additional rotations described with sub-section 5.4.3 creates under column of elements
nonnull under this coefficient. These elements are not awkward because they all are null with the tolerance
chosen.
6 Bibliography
[1]
P.G. CIARLET: “Introduction to the matric numerical analysis and optimization” _ MASSON
(1985).
[2]
G.H. GOLUB, C. REINSCH “Singular been worth decomposition and least public gardens solutions” in
“Handbook for Automatic Computation - Linear Algebra, vol. 2” J.H. WILKHINSON,
C. REINSH Editors _ SPINGER VERLAG (1971).
[3]
P. LASCAUX, R. THEODOR: “Matric numerical Analysis applied to the art of the engineer”,
volumes 1 and 2 _ MASSON (1986).