Regularization of singular least squares problems

Regularization of singular least squares problems

P. Carrette

Department of Electrical Engineering Linkping University, S-581 83 Linkping, Sweden

WWW:

Email:

March 11, 1998

REGLERTEKNIK

AUTOMATIC CONTROL LINKÖPING

Report no.: LiTH-ISY-R-2019 Submitted to BIT journal

Technical reports from the Automatic Control group in Linkping are available by anony-

mous ftp at the address

. This report is contained in the com-

pressed postscript le

.

Regularization of singular least squares problems

P. Carrette

Department of Electrical Engineering, Linkoping University S-58183 Linkoping, Sweden

Abstract

In this note, we analyze the inuence of the regularization procedure applied to singular LS square problems.

It appears that, due to nite numerical accuracy within the computer calculations, the regularization parameter has to belong to a particular range of values in order to have the regularized solution close to that associated to the singular LS problem. Surprisingly enough, this range essentially depends on the square root of the computer precision while the deciency (or singularity) of the regularized LS problem is governed by this precision.

The analysis is based on matrix perturbation theory for which the paper 12] is an utmost reference.

Keywords

: Matrix perturbation, Tikhonov regularization, Singular value decompo- sition.

1 Introduction

In this contribution, we present results concerning the use of Tikhonov regularization (see 13, 10] and references therein) while solving singular least squares (LS) problems, i.e.

x

0

= min

x

kAx;bk

2

subject to min

(1)

for which the matrix

(for

) is (column) rank decient. The corresponding regularized LS problem is as follows

~

x



= min

x

kAx;bk 2

2

+

(2)

for some

value.

Here, we intend to investigate the inuence of the regularization parameter

upon the deviation between the regularized solution ~

and the vector

(solution to problem (1)), i.e. ~

. The reason for this study is that we want to nd a good approximation of

x

0

without solving the singular LS problem itself (by, e.g., performing the pseudo-inverse

1

10⁻¹² 10⁻¹⁰ 10⁻⁸ 10⁻⁶ 10⁻⁴ 10⁻² 10⁰ 10⁻⁸

10⁻⁶ 10⁻⁴ 10⁻² 10⁰ 10² 10⁴ 10⁶ 10⁸

Figure 1: Quantities

~

(

) and

k^x

~

(

) as functions of .

2.2 2.25 2.3 2.35 2.4 2.45 2.5 2.55 2.6 2.65 2.7

10¹ 10²

kA^x~^;bk2

k~xk2

Figure 2:

-curve example with the three points appearing in Figure 1.

11] of the matrix

, i.e.

=

).

As it appears in the literature (see 1, 5, 7, 4, 14] and the Matlab package presented in 8]), Tikhonov regularization gives us a grip for such an approximation while solving the ordinary (full column rank) LS problem (2). The available results about Tikhonov regularization generally consider the links between the two contributions to the cost of the regularized problem, i.e.

and

, only. Indeed, roughly speaking, the solution to the original problem (1) can be obtained in making these two quantities simultaneously small for appropriate value of the regularization parameter

. By this, we have in mind the reasoning that leads to the selection of \best"

values by inspecting the

-curve that is the representation of these two quantities at regularized solutions, i.e.

~

as a function of

~

for dierent values of

, (see 7] for more details).

Unfortunately, the \corner" in the

-curve does not provide

values that are robust with respect to the 2-norm of the deviation we are interested in, i.e.

~

. Let us illustrate this fact by a simple example:

is a 50

5 matrix with rank 4 and unit nonzero singular values while the 50 elements of the column

are samples of an uniformly distributed (in 0

1]) random variable. In Figure 1, we have presented the 2-norm of the deviation between the regularized solution ~

and the original solution

as well as the 2-norm of the regularized solution ~

as functions of the parameter

. The graph of the primer deviation 2-norm can be divided into two parts in agreement with the decreasing and increasing behaviors of this error with

. Three points have been enlightened on the straight line curve. In Figure 2, we have displayed the corresponding

-curve, i.e.

(

~

~

) for dierent values of

. The preceding three points are concentrated in the lower-left corner of this curve so that it does not succeed in giving a high preference to the \o" point that is associated to the \best" approximation of the reference solution

x

0

, i.e. for

5

6 10

.

In order to overcome the poor capability of the

-curve (as well as of the usual studies of Tikhonov regularization) to end up with a

value corresponding to a solution ~

close

2

to the reference vector

, we here provide the analysis of the deviation between these two solutions, i.e. ~

, on the basis of matrix perturbation theory. Therefore, we make an intensive use of the paper of Stewart 12] that presents a complete analysis of the perturbation of the singular value decomposition (SVD) of matrices. Note that in 5], Hansen provides expressions for the upperbound on the relative deviation between a perturbed and the unperturbed solution, denoted

, of a LS problem similar to (2), i.e.

k^x

~

.

As a by-product of our analysis, intervals over the regularization parameter values are found for any given admissible accuracy imposed on the regularized solution ~

, i.e.

] implying that

~

.

The paper is organized as follows. In Section 2, we introduce the notations we will be dealing with in the sequel and we discuss on how numerical errors should be taken into account in order to make the analysis of the regularized solution deviation explain its behavior (as seen in Figure 1). In Section 3, we decompose the deviation between the regularized solution ~

and the original

into two components, i.e. either in or out of the kernel of the problem matrix

. In Section 4, we gives expressions for the singular value decomposition (SVD) of the perturbed version of the matrix

that enters in the resolution of the regularized LS problem (2). This is achieved on the basis of results presented in 12].

In Section 5, we end up with a closed form expression of the 2-norm of the regularized solution error. Its intrinsic characteristics with respect to the regularization parameter

are commented in some details. Finally, simulation examples are presented in Section 6.

They completely agree with the results derived in the preceding sections.

2 Notations and numerical error discussion

First, let us introduce notations for the two LS problems, i.e. minimization (1) and (2), respectively.

The matrix

has a column rank identical to

while

=

+

where

ker

originates from the inconsistency of the LS problem (1). We also dene the matrix

A 2R

(N+n) n

as

A

=



A

0



and

=



b

0



The SVD of

is written as

A

= 

]



S

1

0 0 0



V T

1

V T

2



=

(3)

where

= diag(

) with

, the

-th (in decreasing order) nonzero singular value of

while the matrices 

] and 

] are orthogonal of dimension (

+

) and

, respectively. Then, the orthogonal projector onto the range of

is

=

while its orthogonal counterpart, i.e. the orthogonal projector onto the kernel of

, is denoted

3

K

?

=

(=

).

With these notations in mind, we can write the solution to the original LS problem as

x

0

=

where

is the pseudo-inverse of

, i.e.

=

(uniquely dened). It is worth mentioning that

=

so that

all lies in the range of

.

Remark 1 The \

" operator designed in Matlab for solving LS problems, i.e.

, leads to solutions that generally contains a component in the kernel space of

(when it is non-trivial), i.e.

x

LS

=

+

with nonzero

. The reason for this is that the procedure uses a

decompo- sition of the matrix

and, in case

has its (

) last rows that contain negligible elements regarding the numerical accuracy of the computer, the (

) last elements of the associated LS solution are xed to zero. In fact, this corresponds to choosing

so that

exactly compensates the (

) last elements of the reference solution

(after permutations if necessary in the

decomposition process).

Now, let us turn to the solution of the regularized LS problem presented in (2). This problem can be expressed in terms of a particular perturbation of the matrix

. Obviously, the second term in the brackets of (2) makes the bottom part of the corresponding matrix

A

, denoted

, become an identity matrix scaled by the regularization parameter

, i.e.

A



=



A

I



The corresponding solution is

x



=

(or

)

It is unique because

is full column rank for

0. From a numerical point of view, we must ask for

(e.g.

10

) where

denotes the accuracy of the computer.

Unfortunately, when we simulate the deviation between this \theoretically" regularized solution

and the solution

, it is impossible to explain the deviation generally associ- ated with the \practically" regularized solution (see the straight line in Figure 8 compared to that in Figure 1). This means that our assumption concerning the inuence of the reg- ularization upon the matrix

(leading to

) is not correct. This can be viewed as a bad model of the regularization eect because it cannot reveal its intrinsic characteristics. It actually appears that additional perturbations must be considered: namely, the numerical errors associated to the computation of the regularized solution. More precisely, we deal with the following regularized matrix

~

A



=

+

(4)

4

where

denotes the numerical accuracy, i.e.

= 2

2 10

in Matlab 5.1, and

R

( N+n) n

stands for the (normalized) numerical error matrix (whose structure is detailed below). The SVD of ~

is denoted

~

A



= ~

~

~

]

2

4

~

S

1

0 0 ~

0 0

3

5

 ^V

~

~

V T

2



= ~

~

~

+ ~

~

~

where ~

is a diagonal containing the

largest singular values of ~

(in increasing order!) and ~

is a diagonal containing its (

) last nonzero singular values. The matrices ~

~

~

] and ~

~

] are orthogonal of dimension (

+

) and

, respectively.

The solution of the corresponding LS problem is written as

~

x



= ~

(or ~

) (5)

It is unique in case ~

is full column rank, i.e.

for instance. The purpose of the paper is then to analyze in details the deviation between this solution ~

and the reference solution

dened for the original LS problem.

Finally, let us give a reasonable structure for the numerical error matrix

. Therefore, we consider the left singular subspace associated to

and we show how the associated singular values may induce a particular scaling of the elements of this matrix. We can write

~

A T

 U

1

=

+

(6)

The

-th column of this matrix has a 2-norm that is approximately identical to

. This means that the related numerical perturbation must take this scale into account, i.e.

(

)

+

(

)

=

(

)

+

(

)

for an appropriate matrix

whose denition is made regardless of the singular values of

. This leads to

=

. For what concerns the

counterpart, it can only be said that the largest elements within

will induce large contributions to this matrix product. Hence, we globally propose that the numerical error matrix

can be written as

E

=

+

(7) for which we point out that the matrices

and

have normalized elements (indepen- dently of the

's and of

).

3 Deviation of the regularized solution

The deviation between the regularized solution ~

and the reference solution

can be decomposed into two parts, i.e.

~

x



;x

0

=

~

+

~

(8)

5

The rst term in the right hand side belongs to the range of the transpose of the matrix

while the second completely lies in its kernel (see denition of the orthogonal projector

and

, respectively). Because of these orthogonal projectors, each of these contributions is orthogonal to the other, i.e.

k^x

~

=

~

+

~

In order to analyze these contributions, let us develop the links that exist between the two solution vectors.

From its denition in equation (5), the regularized solution ~

can be written as

~

x



= ~

~

]

 ^S

~

0 0 ~



;2

 ^V

~

~

V T

2

 ^A

~



=

~

~

~

+ ~

~

~



;

V

1 S

2

1 V

T

1 x

0

+

where we have used the fact that

~

A T





=

+

=

+

with

=

. With the structure we have introduced for the numerical error matrix

, we also have that

E T



= (

+

)

=

+

(

)

where it is worth noting that

(identical to

) is the value of the original LS cost at

.

While coming back to the two contributions to the deviation of this regularized solution, we can write

K

k^x

~

=

(

~

)~

(~

)

+

V

1

(

~

)~

(~

)

+

(~

) (

)

+

V

1

(

~

)~

(~

) +

(~

)]

+

(~

) (

)

(9) and

K

?^x

~

=

(

~

)~

(~

) +

(~

)]

+

(~

) (

)

+

V

2

(

~

)~

(~

) +

(~

)]

+

(~

) (

)

(10) Let us give an interpretation for these two expressions. Because of the regularization of the singular LS problem, the right singular vectors of

rotate a little bit, i.e. leading to ~

~

], and its two singular value subsets (i.e.

for

= 1

as well as the remaining

6

~

K k

~

K

?

K k K

?

~

x



x

0

Figure 3: Schematic representation of the compo- nents of the regularized solution ~

.

zero singular values) are also slightly altered giving rise to ~

and ~

. The regularized solution ~

is naturally expressed in terms of these perturbed singular pairs, i.e. (~

~

) and (~

~

). In other words, two components are found for it according to the related sub- spaces for which the orthogonal projectors are ~

= ~

~

and ~

= ~

~

, respectively.

Hence, the (above) expressions for the contributions to the deviation of the regularized solution show up the projection of this solution back to the subspaces associated to the original matrix

, i.e. by use of the orthogonal projector

and

. In Figure 3, we have schematically drawn this back projection for a 2-dimensional case.

It is also worth noticing that the components in the perturbed subspaces are inuenced by the inverse of the square of the corresponding singular values, i.e. ~

and ~

, respec- tively. As the latter will be seen to behave similarly to

, the corresponding component will show an extreme sensitivity to this regularization parameter, i.e.

1

.

In order to analyze these expressions, we must evaluate the role of the regularization parameter

and of the numerical error matrix

on the following quantities

V T

1 ^V

~

~

~

and

~

as well as ~

and ~

. A simple way to achieve this goal is to use results concerning the SVD of perturbed matrices.

4 SVD of perturbed matrices

In 12], Stewart shows that the right singular vectors of a perturbed matrix, ~

say, can be expressed in terms of those of the original matrix,

say, as follows

~

~

] = 

]



I ;P T

P I



(

+

)

0 0 (

+

)



(11) while

is a matrix satisfying the equation system

(

Q

(

+

)

(

+

)

= (

+

)

P

(

+

)

(

+

)

=

(

+

)

(12)

7

for

and 

=

0

]

while

=

.

Note that 

=

and 

= 0 because of (0

]

)(0

]

)

=

together with the orthogonality of the original right singular vectors, i.e. 

].

From expression (11), we immediately have that

V T

1 ^V

~

= (

+

)

~

=

(

+

)

V T

2

~

V

1

=

(

+

)

~

= (

+

)

(13) He also states that the perturbed singular values belong to disjoint sets so that

i

(~

) =

(

+

)

(

+

(

+

))(

+

)

i

(~

) =

(

+

)

(

+

(

))(

+

)

(14) where

(

) denotes the

-th singular value of the matrix

(in increasing order).

After straightforward derivations, we end up with the following result.

Proposition 1 Under the assumption that

, we have that the solution for the

P

and

matrices in the equation system (12) respectively satisfy

P  V T

2 X

T

1

and

where the symbol \

" should be understood in the spectral sense, i.e.

is equivalent to

1. Furthermore, we have that

V T

1

~

V

1

I V T

2

~

V

2

I

and

~

;

~

]

as well as

(~

)

1 +

2

and

(~

)

for appropriate

.

Proof: In case of small perturbations, the last quadratic terms in equations (12) are generally not considered. Thus, we only have to solve

(

Q

(

+

)

(

+

)

= (

+

)

P

(

+

)

(

+

)

=

As the numerical accuracy

is negligible compared to the diagonal elements in

, we get

(

QS

1

;

(

+

)

= (

+

)

P

=

+ (

+

)

Then, we can develop the rst equation as follows

QS 2

1

;

(

+

)(

+

)

=

(

+

) +

(

+

) leading to

QS 2

1

;

2

(

+

)

8

while considering

. Hence, for

, we end up with

(

+

)



as well as

P E T

12 S

;1

1

+ (

+

)(

+

)

where we used the fact that 

= 0. So, from the structure of the numerical error matrix

in expression (7), we get that

P 

(

+

]

)

=

which means that the approximation of

is expressed regardless of the singular values of the original matrix

.

Finally, we are able to evaluate the expressions (13) and (14), respectively.

Therefore, we rst consider approximations of the square root matrices as follows (

+

)

+

2

+

2

(

+

)

+

2

+

2 as 

=

, and

(

+

)

+

2

+

(

)

2

(

+

)

+

2

+

(

)

2

Moreover, we have

(

+

)

(

+

)

(

+

2)

and (

+

)

(

+

)

(

+

2)

=

because 

= 0. Hence, from the denition (14), these expressions lead to approxi- mations of the diagonal elements of the perturbed singular value matrices, i.e.

i

(~

)

1 +

2

and

(~

)

as the diagonal elements of ~

are increasingly ordered and 

=

.

From this proposition, it is immediately seen that the right singular vectors are robust regarding the perturbation related to the regularization parameter

. In fact, this is a well known result concerning the robustness of eigenvectors of symmetric matrices additively perturbed by a scaled identity matrix, e.g.

for instance.

9