Core problém v lineárn´ı aproximaˇcn´ı úloze Ax ≈ b s jednou pravou stranou

(1)

Core probl´ em v line´ arn´ı aproximaˇ cn´ı ´ uloze Ax ≈ b

s jednou pravou stranou

Bakal´ aˇ rsk´ a pr´ ace

Studijn´ı program: B1101 – Matematika Studijn´ı obor: 1101R016 – Matematika Autor pr´ace: Filip J´agr

Vedouc´ı pr´ace: Martin Pleˇsinger

(2)

The core problem within a linear approximation problem Ax ≈ b,

with the single right-hand side

Bachelor thesis

Study programme: B1101 – Mathematics Study branch: 1101R016 – Mathematics

Author: Filip J´agr

Supervisor: Martin Pleˇsinger

(3)

(4)

(5)

Prohl´ aˇ sen´ı

Byl jsem seznámen s t´ım, ˇze na mou bakaláˇrskou práci se plnˇe vztahuje zákon ˇc. 121/2000 Sb., o právu autorském, zejména § 60 – ˇskoln´ı d´ılo.

Beru na vˇedom´ı, ˇze Technická univerzita v Liberci (TUL) neza- sahuje do mých autorských práv uˇzit´ım mé bakaláˇrské práce pro vnitˇrn´ı potˇrebu TUL.

Uˇziji-li bakaláˇrskou práci nebo poskytnu-li licenci k jej´ımu vyuˇzit´ı, jsem si vˇedom povinnosti informovat o této skuteˇcnosti TUL;

v tomto pˇr´ıpadˇe má TUL právo ode mne poˇzadovat úhradu náklad˚u, které vynaloˇzila na vytvoˇren´ı d´ıla, aˇz do jejich skuteˇcné výˇse.

Bakaláˇrskou práci jsem vypracoval samostatnˇe s pouˇzit´ım uvedené literatury a na základˇe konzultac´ı s vedouc´ım mé bakaláˇrské práce a konzultantem.

Souˇcasnˇe ˇcestnˇe prohlaˇsuji, ˇze tiˇstˇen´a verze pr´ace se shoduje s elek- tronickou verz´ı, vloˇzenou do IS STAG.

Datum:

Podpis:

(6)

Anotace

Tato bakaláˇrská práce se zabývá problematikou ˇreˇsen´ı lineárn´ıch aproximaˇcn´ıch úloh ve smyslu tzv. úplných nejmenˇs´ıch ˇctverc˚u (TLS, z anglického total least squares).

V úvodu se seznám´ıme s velmi d˚uleˇzitým nástrojem, kterým je singulárn´ı rozklad (SVD, z anglického singular value decomposition). Dále v práci pop´ıˇseme moˇzné pˇr´ıstupy k ˇreˇsen´ı lineárn´ıch aproximaˇcn´ıch úloh. Nejprve struˇcnˇe zopakujeme bˇeˇznˇe známý pˇr´ıstup pˇreformulován´ı úlohy na tzv. (klasický) problém nejmenˇs´ıch ˇctverc˚u, neboli lineárn´ı regresi. V kontrastu k tomu zavedeme úplný problém nejmenˇs´ıch ˇctverc˚u, neboli tzv. ortogonáln´ı regresi. Dokáˇzeme, ˇze ne vˇzdy má lineárn´ı apro- ximaˇcn´ı úloha ˇreˇsen´ı ve smyslu TLS. Tento d˚ukaz provedeme pomoc´ı ortogonáln´ı transformace p˚uvodn´ı úlohy, tj. zmˇeny báze souˇradného systému, ve kterém je úloha zformulovaná, na blokovˇe diagonáln´ı tvar. Tato transformace povede pˇr´ımo k definici tzv. core problému. Zm´ın´ıme, ˇze core problém lze zapsat ve dvou typických tvarech, pomoc´ı SVD tvaru a pomoc´ı bidiagonáln´ıho tvaru.

V závˇereˇcné ˇcásti pop´ıˇseme software vytvoˇrený v prostˇred´ı Matlab, který dokáˇze generovat pseudonáhodné lineárn´ı aproximaˇcn´ı úlohy pˇredepsaných rozmˇer˚u obsahuj´ıc´ı core problém pˇredepsaných vlastnost´ı. Smyslem tohoto softwaru je vytvoˇren´ı databáze úloh pro statistické testy. To vˇsak jiˇz nen´ı souˇcást´ı této práce.

Kl´ıˇ cov´ a slova:

singulárn´ı rozklad (SVD); (klasický) problém nejmenˇs´ıch ˇctverc˚u (LS); úplný prob- lém nejmenˇs´ıch ˇctverc˚u (TLS); ortogonáln´ı regrese; ortogonáln´ı transformace; core problém; SVD tvar core problému; bidiagonáln´ı tvar core problému

(7)

Abstract

This bachelor thesis deals with solving linear approximation problems in the sense of the so-called total least squares (TLS).

In the beginning we introduce one of the most important tools, the singular value decomposition (SVD). Then we briefly describe possible approaches to solving such problem. First approach is the commonly known (ordinary) least squares formulation (or problem) also know as linear regression. We are particularly interested in the second approach, the total least squares formulation (or problem) also know as orthogonal regression. We show that the linear approximation problem may not have a solution in the sense TLS. This can be show using orthogonal transformations of the original problem, i.e., by changing coordinate systems in which the problem is formulated, so that the transformed problem has a block diagonal form. This transformation leads directly to the definition of the so-called core problem. We mention that the core problem can be written in two typical forms, in the SVD form and in the bidiagonal form.

In the final section, we describe a software tool that we developed in Matlab, which can generate pseudo-random linear approximation problems with given di- mensions and containing a core problem with prescribed properties. The aim of this tool is to assemble a database of such problems, which will be used for futher statistcal testing. This, however, is not part of this thesis.

Key words:

singular value decomposition (SVD); (ordinary) least squares problem (LS); total least squares problem (TLS); orthogonal regression; orthogonal transformation; core problem; SVD form of a core problem; bidiagonal form of a core problem

(8)

Podˇ ekov´ an´ı

Chtˇel bych podˇekovat Martinu Pleˇsingerovi za veden´ı mé bakaláˇrské práce, cenné rady, trpˇelivost a ochotu, kterou mi vˇenoval v pr˚ubˇehu zpracován´ı této práce.

(9)

Obsah

Anotace 5

Abstract 6

Seznam obr´azk˚u 10

Znaˇcen´ı a z´akladn´ı pojmy 11

Uvod´ 13

1 Singul´arn´ı rozklad 16

1.1 Spektr´aln´ı rozklad symetrick´e matice . . . 16

1.2 Symetrick´a matice jako line´arn´ı zobrazen´ı . . . 18

1.3 Obecn´a matice jako line´arn´ı zobrazen´ı . . . 20

1.4 Singul´arn´ı rozklad . . . 21

1.5 Aproximace matice matic´ı niˇzˇs´ı hodnosti . . . 24

2 Upln´´ y probl´em nejmenˇs´ıch ˇctverc˚u 25 2.1 Motivace . . . 25

2.2 Zaveden´ı a odvozen´ı ˇreˇsen´ı pomoc´ı SVD . . . 26

2.3 Existence a jednoznaˇcnost ˇreˇsen´ı. . . 28

3 Core problém 31 3.1 Ortogonáln´ı transformace úlohy . . . 31

3.2 Dovˇetek k vˇetˇe 8 . . . 32

3.3 Zaveden´ı core probl´emu. . . 33

3.4 Tvary core probl´emu . . . 34

3.4.1 Core probl´em v SVD tvaru . . . 34

3.4.2 Bidiagon´aln´ı tvar core probl´emu . . . 35

4 Automatické generován´ı aproximaˇcn´ıch úloh s core problémem s pˇre- dem známými parametry 36 4.1 Program cpSVD. Core problém v SVD tvaru . . . 36

4.2 Program cpBDG. Core problém v bidiagonáln´ım tvaru . . . 41 4.3 Program approxpb. Generátor aproximaˇcn´ıch úloh s core problémem 44

(10)

4.4 Databáze náhodných aproximaˇcn´ıch úloh s core problémem stejných vlastnost´ı . . . 46

Z´avˇer 47

Reference 48

(11)

Seznam obr´ azk˚ u

2.1 Klasický problém nejmenˇs´ıch ˇctverc˚u (lineárn´ı regrese) . . . 26 2.2 Upln´´ y problém nejmenˇs´ıch ˇctverc˚u (ortogonáln´ı regrese). . . 27

(12)

Znaˇ cen´ı a z´ akladn´ı pojmy

Matice a vektory

A ∈ R^n×m (C^n×m) reálná (komplexn´ı) matice m-krát-n s prvky a_i,j A^T transponovaná matice k matici A

A komplexnˇe sdruˇzená matice k matici A A^H = A^T hermitovsky sdruˇzená matice k matici A rank(A) hodnost matice definovaná jako poˇcet lineárnˇe

nezávislých ˇrádk˚u, respektive sloupc˚u matice A Q^T = Q⁻¹ ∈ R^n×n ortogonáln´ı matice

Q^H = Q⁻¹ ∈ C^n×n unit´arn´ı matice

x ∈ Rⁿ (Cⁿ) reálný (komplexn´ı) vektor s délkou n a prvky x_j kxk eukleidovská norma vektoru, kxk = (P

j|x_j|²)^1/2 hx, yi standardn´ı skal´arn´ı souˇcin, hx, yi = y^Hx =P

jx_jy_j x ⊥ y vektory x a y jsou ortogon´aln´ı, tj. hx, yi = 0

kAk₂ spektr´aln´ı norma matice, kAk₂ = max_kxk=1kAxk kAk_F Frobeniova norma matice, kAk_F = (P

i

P

j|a_i,j|²)^1/2 Σ(A) mnoˇzina vˇsech singul´arn´ıch ˇc´ısel matice A vˇcetnˇe

n´asobnost´ı N (A) j´adro matice A

R(A) obor hodnot matice A

dim(V) dimenze lineárn´ıho vektorového prostoru V V^⊥ ortogonáln´ı doplnˇek podprostoru V

V ⊥ W podprostory V a W jsou ortogon´aln´ı, tj. hx, yi = 0, ∀x ∈ V a ∀y ∈ W V ⊕ W je direktn´ı souˇcet podprostor˚u V a W span(q₁, . . . , q_m) prostor generovan´y vektory q₁, . . . , q_m min M, M ⊂ R minimum z mnoˇziny M

|M| poˇcet prvk˚u koneˇcn´e mnoˇziny M

(13)

Terminologie aproximaˇ c´ıch ´ uloh

Ax ≈ b aproximaˇcn´ı ´uloha s matic´ı A ∈ R^n×m a vektorem prav´e strany b ∈ Rⁿ, b /∈ R(A)

A₁₁x₁ ≈ b₁ core problém v úloze Ax ≈ b, A₁₁ ∈ R^n×m, b₁ ∈ Rⁿ A systémová matice (matice modelu)

b pravá strana (vektor pozorován´ı) [b, A] rozˇs´ıˇrená matice (matice dat)

Pouˇ zit´ e zkratky

SVD singul´arn´ı rozklad matice, A = U ΣV^H (singular value decomposition)

TLS úplný problém nejmenˇs´ıch ˇctverc˚u (total least squares problem)

(14)

Uvod ´

Ustˇredn´ım bodem t´´ eto pr´ace bude aproximaˇcn´ı probl´em Ax ≈ b,

kde A ∈ R^n×m je daná matice, reprezentuj´ıc´ı napˇr. nˇejaký model a pravá strana b ∈ Rⁿ je daný vektor pozorován´ı, který m˚uˇze být pˇrirozenˇe zat´ıˇzený chybami, napˇr. chybami mˇeˇren´ı. Obecnˇe tedy vektor b nemus´ı leˇzet v oboru hodnot R(A) matice A. Proto budeme dále striktnˇe pˇredpokládat

b 6∈ R(A),

v opaˇcném pˇr´ıpadˇe by byl aproximaˇcn´ı problém de-facto soustavou rovnic Ax = b a mˇel by ˇreˇsen´ı v klasickém slova smyslu. Protoˇze ale problém Ax ≈ b, b 6∈ R(A) ˇreˇsen´ı nemá (neexistuje vektor x takový, aby nastala rovnost), mus´ıme pˇristoupit k nˇejakému náhradn´ımu postupu. Typicky hledáme vektor x takový, abychom mini- malizovali nˇejaký pˇredepsaný funkcionál. ˇCasto se v takovém pˇr´ıpadˇe pouˇz´ıvá nˇejaká varianta metody nejmenˇs´ıch ˇctverc˚u. V této práci se budeme zabývat tzv. úplným problémem nejmenˇs´ıch ˇctverc˚u (TLS, z anglického total least squares), který hledá nejmenˇs´ı matici E a vektor r, tj.

min k[r, E]k_F tak, ˇze (b + r) ∈ R(A + E),

kde k · k_F znaˇc´ı tzv. Frobeniovu normu matice, tj. odmocninu ze souˇctu kvadrát˚u (absolutn´ıch hodnot) vˇsech prvk˚u matice. Pokud taková minimáln´ı E a r existuj´ı, pak libovolný vektor x splˇnuj´ıc´ı (A + E)x = (b + r) nazýváme TLS ˇreˇsen´ı p˚uvodn´ı

´

ulohy Ax ≈ b. Tato úloha (na rozd´ıl od klasického problému nejmenˇs´ıch ˇctverc˚u, kde hledáme pouze minimáln´ı opravu r pravé strany b) obecnˇe nemá ˇreˇsen´ı pro dané A a b, jak jiˇz v roce 1980 ukázali G. H. Golub a C. F. Van Loan v ˇclánku [8].

Radu praktick´ˇ ych aplikac´ı úplného problému nejmenˇs´ıch ˇctverc˚u nalezneme napˇr.

v knih´ach [19], [20].

Dalˇs´ı významný posun v analýze TLS problému byl publikován v roce 1991 v knize [18] J. Vandewalle a S. Van Huffel (vycházej´ıc´ı z dizertaˇcn´ı práce S. Van Huffel [15]). Zde se mimo jiné zavád´ı klasifikace problém˚u Ax ≈ b na tzv. gene- rické (generic) a negenerické (nongeneric) problémy, pˇriˇcemˇz prvn´ı z nich vˇzdy maj´ı TLS ˇreˇsen´ı (ne nutnˇe jednoznaˇcné), zat´ımco negenerické problémy TLS ˇreˇsen´ı nemaj´ı. Term´ın

”generický/negenerický problém“ ovˇsem nen´ı nejˇst’astnˇeji zvolen. Po- jem ”generic“ lze z angliˇctiny pˇreloˇzit jako

”bˇeˇzn´y“ a jeho antonymum

”nongeneric“

(15)

pak napˇr. jako

”vzácný“. Tato terminologie byla pravdˇepodobnˇe motivována t´ım, ˇze se bˇeˇzných výpoˇcetn´ıch úlohách skuteˇcnˇe vyskytovaly zpravidla jen generické problémy.

Nejjednoduˇsˇs´ı pˇr´ıpad TLS problému, který je v literatuˇre také nejˇcastˇeji zmiˇno- ván, nav´ıc pracuje s matic´ı A plné sloupcové hodnosti, tj. rank(A) = m. Spojen´ı tˇechto dvou jev˚u nás m˚uˇze snadno svést k následuj´ıc´ı myˇslence:

”Má-li matice A lineárnˇe nezávislé sloupce, pak má problém Ax ≈ b TLS ˇreˇsen´ı,“ která vˇsak nen´ı pravdivá. S. Van Huffel v ˇclánku [17] na str. 547 opatrnˇe p´ıˇse:

”These conditions (podm´ınky postaˇcuj´ıc´ı pro existenci TLS ˇreˇsen´ı) are generically satisfied provided A is of full rank and the set Ax ≈ b is not too conflicting.¹“ V dalˇs´ıch pracech m˚uˇzeme nalézt i ménˇe opatrné formulace, napˇr.:

”If we suppose that [b, A] has full column rank, this condition (podm´ınka postaˇcuj´ıc´ı pro existenci TLS ˇreˇsen´ı) is generally satisfied,“ viz [5, str. 54]; nebo dokonce pˇr´ımo nepravdiv´e, napˇr.:

”In our case, when the matrix A has full rank and n > m, the solution of the TLS problem always exists,“ viz [6, str. 37].

Nepravdivost výˇse zm´ınˇeného tvrzen´ı, tedy fakt, ˇze existuj´ı problémy Ax ≈ b, b 6∈ R(A), rank(A) = m, které nemaj´ı TLS ˇreˇsen´ı, lze nejsnáze ukázat pomoc´ı tzv. core problému zavedeného C. C. Paigem a Z. Strakoˇsem v ˇclánku [12]. Uˇzit´ım core problému lze dokonce takové problémy snadno konstruovat. C´ılem této práce je vytvoˇrit software v prostˇred´ı Matlab, který bude generovat pseudonáhodné apro- ximaˇcn´ı problémy Ax ≈ b (tedy matice A a vektory b), pˇredepsaných rozmˇer˚u, obsahuj´ıc´ı core problém pˇredepsaných rozmˇer˚u a maj´ıc´ı dalˇs´ı pˇredepsané vlastnosti.

T´ım budeme schopni generovat velké mnoˇzstv´ı aproximaˇcn´ıch úloh s maticemi plné sloupcové hodnosti a nemaj´ıc´ı TLS ˇreˇsen´ı.

Dalˇs´ım krokem, který by mˇel navazovat, ale nen´ı jiˇz souˇcást´ı této práce, bude vy- tvoˇrit nˇekolik datových soubor˚u s problémy Ax ≈ b (kaˇzdý soubor charakterizovaný parametry, podle kterých byl vytvoˇren) dostateˇcnˇe velkých pro následuj´ıc´ı statistic- kou analýzu. Konkrétnˇe by bylo zaj´ımavé nauˇcit se generovat (tj. odladit parametry) problémy Ax ≈ b, jejichˇz prvky (tj. prvky jejich matic A a pravých stran b) se budou chovat jako nezávislé, identicky distribuované náhodné promˇenné, matice A by mˇely lineárnˇe nezávislé sloupce, a problémy by pˇritom nemˇely TLS ˇreˇsen´ı. Zaj´ımavý by byl i negativn´ı výsledek, tedy pokud by se takové problémy generovat nepodaˇrilo.

Naznaˇcovalo by to, ˇze existenci core problému uvnitˇr aproximaˇcn´ıho problému by mohlo být moˇzné detekovat statistickými metodami.

Struktura této práce je následuj´ıc´ı. Po struˇcném úvodu, v kapitole 1 zavedeme tzv. singulárn´ı rozklad matice, jehoˇz zaveden´ı motivujeme spektráln´ım rozkladem.

V z´avˇeru kapitoly uk´aˇzeme, jak lze danou matici aproximovat matic´ı niˇzˇs´ı hodnosti.

V kapitole2se budeme vˇenovat úplnému problému nejmenˇs´ıch ˇctverc˚u (TLS), který motivujeme jednoduchým fyzikáln´ım pˇr´ıkladem na elektrickou vodivost a tzv. kla- sickým problémem nejmenˇs´ıch ˇctverc˚u. V závˇeru kapitoly zformulujeme postaˇcuj´ıc´ı podm´ınku pro existenci ˇreˇsen´ı aproximaˇcn´ı úlohy ve smyslu TLS. V kapitole 3 se

1Cl´ˇ anek [17] je také k dispozici v preprintu [16], kde je citovaná vˇeta na str. 9. V ˇclánku [17] i v jeho preprintu [16] se m´ısto Ax ≈ b pouˇz´ıvá znaˇcen´ı Xβ ≈ y obvyklejˇs´ı ve statistické literatuˇre. Zde je znaˇcen´ı obmˇenˇeno z d˚uvodu konzistence textu. Obdobnˇe je znaˇcen´ı obmˇenˇeno v obou následuj´ıc´ıch vˇetách citovaných z ˇclánk˚u [5] a [6].

(16)

budeme vˇenovat tzv. core problému nástroji, který vˇzdy nalezne ˇreˇsen´ı aproximaˇcn´ı

´

ulohy za pˇredpokladu, ˇze ˇreˇsen´ı existuje. V kapitole se budeme také vˇenovat r˚uzným tvar˚um core problému a sice tzv. SVD tvaru core problému a tzv. Bidiagonáln´ımu tvaru core problému. V kapitole 4 se budeme vˇenovat popisu naprogramovaných program˚u cpSVD a cpBDG, které generuj´ı core problém v SVD respektive v bidia- gonáln´ım tvaru a také programu aproxpb, který bude generovat aproximaˇcn´ı úlohu Ax ≈ b s core problémem.

(17)

1 Singul´ arn´ı rozklad

V následuj´ıc´ı kapitole se budeme vˇenovat singulárn´ımu rozkladu (SVD). Ve výkladu SVD budeme vycházet ze skript [2], respektive z knih [21], [4], [9]. SVD je silným nástrojem pro výpoˇcet mnoha vlastnost´ı matic. Pomoc´ı SVD m˚uˇzeme spoˇc´ıtat napˇr´ıklad hodnost, spektráln´ı nebo Frobeniovu normu matice, pseudoinverzn´ı matici, bázi nulového prostoru a oboru hodnot matice a ˇradu dalˇs´ıch vlastnost´ı. Za variabilitu a s´ılu SVD nicménˇe mus´ıme poˇc´ıtat s vyˇsˇs´ımi výpoˇcetn´ımi nároky.

Odvozen´ı singulárn´ıho rozkladu budeme motivovat spektráln´ım rozkladem sy- metrické matice.

1.1 Spektr´ aln´ı rozklad symetrick´ e matice

Spektráln´ı rozklad symetrické semidefinitn´ı matice popisuje následuj´ıc´ı vˇeta.

Vˇeta 1. Necht’ S ∈ R^n×n je symetrick´a pozitivnˇe semidefinitn´ı matice s hodnost´ı rank(S) = r. Potom existuje ortogon´aln´ı matice Q ∈ R^n×n tak, ˇze

D = Q^TSQ (1.1)

je diagonáln´ı matice s vlastn´ımi ˇc´ısly matice S na diagonále. Ta jsou reálná a seˇrazená tak, ˇze plat´ı

λ₁ ≥ λ₂ ≥ ... ≥ λ_r > λ_r+1 = . . . = λ_n = 0.

D˚ukaz. Schurova vˇeta [2, vˇeta 2.2] ˇr´ıká, ˇze pro kaˇzdou matici S existuje unitárn´ı matice Q ∈ C^n×n, Q = Q^H = Q^T, taková, ˇze

S = QRQ^H, kde

R = Q^HSQ

je komplexn´ı horn´ı trojúheln´ıková matice s vlastn´ımi ˇc´ısly matice S na diagonále seˇrazenými v libovolném poˇrad´ı. D˚ukaz bude veden ve tˇrech kroc´ıch. V prvn´ım kroku ukáˇzeme, ˇze vlastn´ı ˇc´ısla matice S jsou reálná. Dále ukáˇzeme, ˇze vlastn´ı ˇc´ısla jsou nezáporná. Nakonec ukáˇzeme, ˇze matici Q lze volit reálnou.

(18)

Krok 1: Je zˇrejm´e, ˇze

S^T = (QRQ^H)^T = QR^HQ^H.

Protoˇze je matice S podle pˇredpokladu vˇety symetrick´a, tedy S = S^T, plat´ı QRQ^H = QR^HQ^H.

Vyn´asoben´ım obou stran rovnice matic´ı Q^H zleva z´ısk´ame RQ^H = R^HQ^H

dalˇs´ım vyn´asoben´ım obou stran t´eto rovnice matic´ı Q zprava vid´ıme, ˇze plat´ı R = R^H = R^T,

tedy je patrné, ˇze matice R mus´ı být diagonáln´ı a reálná.

Krok 2: Protoˇze matice S je pozitivnˇe semidefinitn´ı, tedy x^TSx ≥ 0

pro libovoln´e x ∈ Rⁿ, pak

q^T_jSq_j = q^T_j(q_jλ_j) = (q_j^Tq_j)λ_j = λ_j ≥ 0 pro j = 1, . . . , n. Protoˇze plat´ı

rank(S) = rank(QDQ^T) a Q je regul´arn´ı, potom

rank(S) = rank(D),

protoˇze matice D je diagon´aln´ı, t.j. D = diag(λ₁, λ₂, ..., λ_n), a rank(D) = r,

pak právˇe n − r vlastn´ıch ˇc´ısel matice A mus´ı být rovno nule. Tedy zbývaj´ıc´ıch r vlastn´ıch ˇc´ısel je kladných.

Krok 3: Vztah S = QDQ^H lze pˇrepsat ve tvaru

SQ = QD, (1.2)

t.j.

Sq_j = q_jλ_j, j = 1, ..., n, (1.3) kde q_j je j-tý sloupec matice Q a λ_j je j-tý diagonáln´ı prvek matice D. Tedy q_j je vlastn´ı vektor matice S odpov´ıdaj´ıc´ı vlastn´ımu ˇc´ıslu λ_j, t.j.

(S − Iλ_j)q_j = 0. (1.4)

Protoˇze matice S − Iλ_j je reálná má homogenn´ı soustava (1.4) také reálné ˇreˇsen´ı.

Tedy vlastn´ı vektor q_j a celou matici Q lze volit re´aln´e.

(19)

1.2 Symetrick´ a matice jako line´ arn´ı zobrazen´ı

Kaˇzdá ˇctvercová matice S ∈ R^n×n pˇredstavuje lineárn´ı zobrazen´ı z prostoru Rⁿzpˇet do prostoru Rⁿ. Toto zobrazen´ı pˇriˇrad´ı obecnému vektoru x ∈ Rⁿ vektor Sx ∈ Rⁿ, tedy

S : x 7−→ y = Sx.

Pˇripomeˇnme, ˇze pro libovolné lineárn´ı zobrazen´ı (a tedy i pro matici) lze zavést pojem obor hodnot zobrazen´ı a jádro zobrazen´ı. Nejprve zavedeme pojem obor hodnot.

Definice 1 (Obor hodnot matice). Necht’ A ∈ R^n×m je obecnˇe obdeln´ıkov´a matice.

Mnoˇzina

R(A) = {y ; y = Ax, x ∈ R^m} ⊆ Rⁿ se naz´yv´a obor hodnot matice A.

N´asleduj´ıc´ı vˇeta upˇresn´ı jakou strukturu m´a obor hodnot matice.

Vˇeta 2. Necht’ A ∈ R^n×m je obecnˇe obdeln´ıková matice. Mnoˇzina R(A) z definice 1 je lineárn´ı vektorový prostor.

D˚ukaz. D˚ukaz je jednoduchý. Uvaˇzujme dva libovolné vektory z oboru hodnot matice, napˇr. y₁, y₂ ∈ R(A), pak existuj´ı vektory x₁, x₂, pro které plat´ı y₁ = Ax₁ a y₂ = Ax₂. Dále zvolme x = αx₁+ βx₂, potom

Ax = A(αx₁+ βx₂) = αAx₁+ βAx₂ = αy₁ + βy₂ = y ∈ R(A).

Nyn´ı zavedeme druh´y z pojm˚u, tedy j´adro matice.

Definice 2 (J´adro matice). Necht’ A ∈ R^n×mje obecnˇe obdeln´ıkov´a matice. Mnoˇzina N (A) = {x ; Ax = 0} ⊆ R^m

se nazývá jádro matice A.

Stejnˇe jako u oboru hodnot, n´asleduj´ıc´ı vˇeta upˇresn´ı strukturu j´adra matice.

Vˇeta 3. Necht’ A ∈ R^n×m je obecnˇe obdeln´ıková matice. Mnoˇzina N (A) z definice 2 je lineárn´ı vektorový prostor.

D˚ukaz. Pro vektory plat´ı x₁, x₂ ∈ R(A). D´ale zvolme x = αx₁+ βx₂, potom Ax = A(αx₁+ βx₂) = αAx₁+ βAx₂ = α0 + β0 = 0.

(20)

Protoˇze matice S je symetrická pozitivnˇe semidefinitn´ı a podle vˇety 1 v´ıme, ˇze má n vlastn´ıch vektor˚u a nav´ıc, ˇze tyto vektory lze volit navzájem ortogonáln´ı. Tedy vlastn´ı vektory q_j, j = 1, . . . , n, matice S tvoˇr´ı ortonormáln´ı bázi Rⁿ. Tedy ze vztah˚u (1.1) respektive (1.2) a (1.3) plyne

S : q1 7−→ λ1 q1, S : q₂ 7−→ λ₂ q₂,

...

S : q_r 7−→ λ_r q_r, S : {q_r+1, . . . , q_n} 7−→ 0 .

(1.5)

Libovolný vektor x ∈ Rⁿ m˚uˇzeme zapsat jako lineárn´ı kombinaci bázových vektor˚u q_j, j = 1, . . . , n následuj´ıc´ım zp˚usobem

x =

n

X

j=1

αjqj. (1.6)

C´ısla αˇ j jsou souˇradnice vektoru x ve smˇeru qj. Pak

S : x 7−→

n

X

j=1

α_j(λ_jq_j).

Z toho plyne

R(S) = span(q1, q2, . . . , qr) = R(S^T),

N (S) = span(q_r+1, . . . , q_n) = N (S^T). (1.7) Protoˇze plat´ı Q^TQ = I vid´ıme, ˇze pro ∀y ∈ R(S) a pro ∀x ∈ N (S) plat´ı, ˇze jejich skal´arn´ı souˇcin bude roven nule, t.j.

hx, yi = 0.

Ze skalárn´ıho souˇcinu plyne, ˇze obor hodnot matice a jádro matice jsou na sebe kolmé, tedy

R(S) ⊥ N (S). (1.8)

Vektor (1.6) lze tak´e zapsat ve tvaru

x = x_R+ x_N, kde x_R a x_N jsou definov´any n´asledovnˇe

x_R =

r

X

j=1

α_jq_j ∈ R(S),

xN =

n

X

j=r+1

αjqj ∈ N (S).

(21)

Pak pro x_R a x_N plyne z (1.8), ˇze

x_R⊥ x_N. Tuto vlastnost zapisujeme

N (S) ⊕ R(S) = Rⁿ, (1.9)

kde operace ⊕ se nazývá direktn´ı souˇcet podprostor˚u. V dalˇs´ı kapitole se pokus´ıme z´ıskat popis analogický k (1.5) a (1.7) pro obecnou obdeln´ıkovou matici.

1.3 Obecn´ a matice jako line´ arn´ı zobrazen´ı

Protoˇze dále budeme pracovat s obecnou obdeln´ıkovou matic´ı je tˇreba vyjasnit následuj´ıc´ı. Ze základn´ıho pˇrehledu pˇrednáˇsek z lineárn´ı algebry v´ıme, ˇze poˇcet lineárnˇe nezávislých ˇrádk˚u se rovná poˇctu lineárnˇe nezávislých sloupc˚u. Tedy pro obecnou obdeln´ıkovou matici A ∈ R^n×m plat´ı

rank(A) = rank(A^T),

pˇriˇcemˇz hodnost matice A je definov´ana jako dimenze prostoru R(A), tedy rank(A) = dim(R(A)).

Zˇrejmˇe tedy plat´ı

dim(R(A)) = rank(A) = rank(A^T) = dim(R(A^T)). (1.10) Plat´ı n´asleduj´ıc´ı vˇeta zobecˇnuj´ıc´ı vztahy (1.8) a (1.9).

Vˇeta 4. Uvaˇzujme obecnou obdeln´ıkovou matici A ∈ R^n×m. Pro kterou plat´ı N (A) ⊕ R(A^T) = R^m, R(A^T) ⊥ N (A), (1.11)

N (A^T) ⊕ R(A) = Rⁿ, R(A) ⊥ N (A^T). (1.12)

D˚ukaz. D˚ukaz pˇrevezmeme z [2, vˇeta 1.25]. Plat´ı R(A^T) ⊕ R(A^T)^⊥ = R^m.

Pro dokázán´ı prvn´ı ˇcásti staˇc´ı ukázat, ˇze R(A^T)^⊥ ⊥ N (A). Zavedeme libovolný vektor x z R(A^T)^⊥ pro který plat´ı následuj´ıc´ı ekvivalence

x ∈ R(A^T)^⊥ ⇐⇒ pro kaˇzd´e y ∈ R(A^T) plat´ı hy, xi = 0

⇐⇒ pro kaˇzd´e z ∈ Rⁿplat´ı A^Tz, x = 0

⇐⇒ pro kaˇzd´e z ∈ Rⁿplat´ı hz, Axi = 0

⇐⇒ Ax = 0

⇐⇒ x ∈ N (A).

Druhá ˇcást tvrzen´ı vyplývá z vyuˇzit´ı prvn´ıho tvrzen´ı na matici A^T.

(22)

Pro matici A a symetrickou pozitivnˇe semidefinitn´ı matice A^TA, AA^T nav´ıc plat´ı n´asleduj´ıc´ı vˇeta.

Vˇeta 5. Pro obecnou obdeln´ıkovou matici A ∈ R^n×m plat´ı N (A^TA) = N (A), R(A^TA) = R(A^T), N (AA^T) = N (A^T), R(AA^T) = R(A).

D˚ukaz. D˚ukaz opˇet pˇrevezmeme z [2, lemma 5.1]. Nejprve dok´aˇzeme prvn´ı rovnost.

Necht’ x ∈ N (A), potom plat´ı Ax = 0 a tak´e A^TAx = 0. Z toho plyne, ˇze x ∈ N (A^TA). Pak tedy

N (A) ⊂ N (A^TA).

Nyn´ı necht’ x ∈ N (A^TA), potom plat´ı AA^TAx = 0 a tak´e

0 = x^TA^TAx =A^TAx, x = hAx, Axi = kAxk². Tedy i Ax = 0 a potom x ∈ N (A). Pak tedy

N (A^TA) ⊂ N (A).

T´ım je dok´az´ana rovnost N (A^TA) = N (A). Rovnost pro obory hodnot, tj. R(A^TA) = R(A^T), dostaneme pˇr´ımo z rovnosti pˇredchoz´ı a ze vztahu (1.11).

Pro d˚ukaz druh´e sady rovnosti staˇc´ı uvaˇzovat transponovanou matici B = A^T, a do prvn´ı sady rovnost´ı dosadit za A i A^T. T´ım ihned dostaneme vztahy N (BB^T) = N (B^T) a R(BB^T) = R(B).

Podle rovnosti (1.10) a podle vˇety 5dost´av´ame

dim(R(AA^T)) = dim(R(A)) = r = dim(R(A^T)) = dim(R(A^TA)), (1.13) kde r je hodnost matice A.

1.4 Singul´ arn´ı rozklad

Následuj´ıc´ı vˇeta o singulárn´ım rozkladu (SVD, z anglického singular value decomposition) bude pro dalˇs´ı výklad kl´ıˇcová.

Vˇeta 6 (O singulárn´ım rozkladu). Necht’ A ∈ R^n×m je obecná obdeln´ıková matice hodnosti r. Potom existuj´ı ortogonáln´ı matice U ∈ R^n×n, V ∈ R^m×m a matice Σ ∈ R^n×m maj´ıc´ı nenulové prvky pouze na diagonále tak, ˇze plat´ı

A = U ΣV^T. (1.14)

Sloupce matic U a V ,

U = [u1, . . . , un], V = [v1, . . . , vm],

(23)

nazýváme levé respektive pravé singulárn´ı vektory. Matici Σ lze zapsat ve tvaru Σ = Σ_r 0

0 0

,

Σ_r =





 σ₁

. ..

σ_r





, pˇriˇcemˇz plat´ı

σ₁ ≥ σ₂ ≥ . . . ≥ σ_r> 0. (1.15) Nenulová ˇc´ısla σj, j = 1, . . . , r nazýváme singulárn´ı ˇc´ısla.

D˚ukaz. Protoˇze A^TA ∈ R^m×m je symetrická pozitivnˇe semidefinitn´ı podle vˇety 1 existuje ortogonáln´ı matice V a diagonáln´ı matice D tak, ˇze plat´ı

D = diag(λ₁, . . . λ_m) = V^T(A^TA)V.

Protoˇze rank(A^TA) = rank(A) = r podle (1.13) a v´ıme, ˇze λ₁ ≥ λ₂ ≥ . . . ≥ λ_r > λ_r+1 = . . . = λ_n = 0.

Tedy plat´ı

A^TAv_j = v_jλ_j pro j = 1, . . . , r a nav´ıc

N (A^TA) = N (A) = span(v_r+1, . . . , v_m), viz (1.5), (1.7). Zˇrejme pro j = 1, . . . , r plat´ı

A(A^TA)v_j = Av_jλ_j, (AA^T)(Av_j) = (Av_j)λ_j,

AA^Tw_j = w_jλ_j, kde

w_j = Av_j

a matice AA^T ∈ R^n×n je symetrick´a pozitivnˇe semidefinitn´ı s hodnost´ı r. Pro skal´arn´ı souˇcin vektor˚u w_j a w_k plat´ı

hw_j, w_ki = w_k^Tw_j = v^T_k(A^TAv_j) = v_k^Tλ_jv_j = λ_jhv_j, v_ki , tedy

hw_j, w_ki =

0, j 6= k λ_j, j = k . Normalizac´ı vektor˚u w_j z´ısk´ame

u_j = w_j

kwjk = 1 pλ_jw_j.

(24)

Protoˇze

dim(R(AA^T)) = r (viz (1.13)), pak

R(AA^T) = span(v₁, . . . , v_r) z vˇety 4vypl´yv´a

dim(N (AA^T)) = n − r, a d´ale

R(AA^T) ⊥ N (AA^T).

Pro libovolnou ortonorm´aln´ı b´azi {u_r+1, . . . , u_n} prostoru N (AA^T) tedy plat´ı, ˇze U = [u₁, . . . , u_r, u_r+1. . . , u_n]

je ortogon´aln´ı matice, t.j. U U^T = U^TU = I_n. Tedy D = U^T(AA^T)U.

Pro j = 1, . . . , r tedy plat´ı

AA^Tu_j = u_jλ_j, A^T(AA^T)u_j = A^Tu_jλ_j, (A^TA)(A^Tu_j) = (A^Tu_j)λ_j,

A^TAwej = wejλj, kde

wej = A^Tuj. Pro skal´arn´ı souˇcin vektor˚u we_j awe_k plat´ı

hwej,weki =

0, i 6= k λ_j, j = k . Normalizac´ı vektor˚uwe_j z´ısk´ame

ev_j = we_j

kwe_jk = 1 pλ_jwe_j. Nyn´ı zbývá ukázat, ˇze vj =evj, zˇrejmˇe

ev_j = 1 pλj

A^Tu_j = 1

λ_jA^Tw_j = 1

λ_jA^TAv_j = 1

λ_jλ_jv_j = v_j. Oznaˇc´ıme-li

σ_j =pλ_j, j = 1, . . . , r pak je d˚ukaz hotov.

Ukázali jsme tedy, ˇze pro libovolnou matici A existuje singulárn´ı rozklad ve tvaru (1.14). V následuj´ıc´ı sekci ukáˇzeme dalˇs´ı moˇzné tvary singulárn´ıho rozkladu a jejich vyuˇzit´ı k aproximaci matice matic´ı niˇzˇs´ı hodnosti.

(25)

1.5 Aproximace matice matic´ı niˇ zˇ s´ı hodnosti

V nˇekterých praktických úlohách je lepˇs´ı pouˇz´ıt tzv. ekonomický tvar singulárn´ıho rozkladu, který spoˇc´ıvá v odstranˇen´ı nˇekterých blok˚u, které pˇri násoben´ı matic U a V^T s matic´ı Σ nemaj´ı vliv na výsledný souˇcin, viz následuj´ıc´ı schéma (pˇrevzato z [2, str. 127]):

A

=

U Σ

Σ_r

U_r V_r^T

0 0 0

V^T

=

U_r Σ_r V_r^T .

Plat´ı tedy

A = U ΣV^T = U_rΣ_rV_r^T. (1.16) Z ekonomick´eho tvaru je vidˇet, ˇze matici A m˚uˇzeme zapsat pomoc´ı souˇctu vnˇejˇs´ıch souˇcin˚u

A =

k

X

j=1

σjujv^T_j

t.j. takzvaný dyadický rozvoj. Následuj´ıc´ı d˚uleˇzitá vˇeta popisuje aproximaci matice matic´ı niˇzˇs´ı hodnosti.

Vˇeta 7 (Schmidt, Eckart, Young, Mirsky). Necht’ A ∈ R^n×m je matice s hodnost´ı r. Pro libovolné pˇrirozené ˇc´ıslo k, pro které plat´ı k < r, matice

A^(k)=

k

X

j=1

σ_ju_jv^T_j

je nejlepˇs´ı aproximac´ı matice A v n´asleduj´ıc´ım smyslu:

min

X∈R^n×m, rank(X)≤kkA − Xk₂ = kA − A^(k)k₂ = σ_k+1 a

min

X∈R^n×m, rank(X)≤kkA − Xk_F = kA − A^(k)k_F = v u u t

r

X

j=k+1

σ_j².

D˚ukaz pro jednotlivé normy lze nalézt v p˚uvodn´ıch ˇclánc´ıch [3] a [11]. D˚ukaz pro libovolnou ortogonálnˇe invariantn´ı maticovou normu lze nalézt v knize [13, str. 208–

209]. D˚ukaz pro spektráln´ı normu je také ve skriptech [2, str. 133, vˇeta 5.12]. Pozna- menejme, ˇze nejlepˇs´ı aproximace z pˇredchoz´ı vˇety nen´ı dána jednoznaˇcnˇe napˇr´ıklad kdyˇz σ_k = σ_k+1.

(26)

2 Upln´ ´ y probl´ em nejmenˇ s´ıch ˇ ctverc˚ u

V této kapitole se jiˇz budeme zabývat ˇreˇsen´ım aproximaˇcn´ı úlohy

Ax ≈ b s podm´ınkou b /∈ R(A) (2.1)

a to pomoc´ı úplného problému nejmenˇs´ıch ˇctverc˚u (TLS, z anglického total least squares). K zaveden´ı TLS budeme vycházet ze skript [2].

2.1 Motivace

Na úvod uved’me jednoduchý pˇr´ıklad na elektrickou vodivost, která je definovaná jako pod´ıl elektrického proudu a napˇet´ı. Z namˇeˇrených hodnot urˇc´ıme jednotlivé vodivosti viz graf na obrázku2.1. Následnˇe metodou nejmenˇs´ıch ˇctverc˚u aproximujeme namˇeˇrená data pˇr´ımkou.

Nejprve zaˇcneme s definic´ı klasick´eho probl´emu nejmenˇs´ıch ˇctverc˚u.

Definice 3 (Klasický problém nejmenˇs´ıch ˇctverc˚u). Necht’ existuj´ı A ∈ R^n×m a b ∈ Rⁿ. Klasickým problémem nejmenˇs´ıch ˇctverc˚u nazveme úlohu nalezen´ı nejmenˇs´ı opravy pravé strany tak, aby opravená úloha byla kompatibiln´ı t.j.

min krk_F tak, ˇze (b + r) ∈ R(A).

Pak existuje vektor x ∈ R^m, který ˇreˇs´ı opravenou úlohu. Problém lze tedy také formulovat následuj´ıc´ım zp˚usobem

min krk_F tak, ˇze Ax = b + r, x ∈ R^m.

Nˇekdy se také m˚uˇzeme setkat, ˇze problém nejmenˇs´ıch ˇctverc˚u je oznaˇcován jako lineárn´ı regrese.

Metoda nejmenˇs´ıch ˇctverc˚u (klasický problém je jej´ı variantou v lineárn´ım pˇr´ıpadˇe;

obecnˇe A(x) nemus´ı být lineárn´ı zobrazen´ı) je pˇripisován C. F. Gaußovi (1777–1855), jeho základy lze nalézt v práci [7]. Ze souˇcasných uˇcebnic zabývaj´ıc´ıch se ˇreˇsen´ım problému nejmenˇs´ıch ˇctverc˚u odkazujeme na [10] a [1]. Z definice a také z pˇr´ıkladu na obrázku2.1 na výpoˇcet vodivosti je patrné, ˇze v modelu je opravována pouze jedna sada namˇeˇrených dat. V pˇr´ıkladu na obrázku2.1jsou to data na ose U (napˇet´ı). Data na ose I (proudy) jsou povaˇzována za pˇresná. Nicménˇe v praxi mohou být obsaˇzeny chyby v obou sadách namˇeˇrených dat. Tomu se budeme vˇenovat v následuj´ıc´ı sekci.

(27)

0 2 4 6 8 10 12 14 16 18 0

2 4 6 8 10 12 14 16 18

U [V]

I [A]

Obrázek 2.1: Body pˇredstavuj´ı namˇeˇrené hodnoty proudu I (svislá osa) procházej´ıc´ıho rezistorem pˇri deseti r˚uzných mˇeˇr´ıc´ıch napˇet´ıch U (vodorovná osa).

Ukolem je urˇ´ cit vodivost G, resp. rezistivitu R souˇcástky (R = G⁻¹). Z Ohmova zákona plat´ı UG = I. Úloha má tedy tvar Ax ≈ b, A ∈ R^n×m, x ∈ R^m, b ∈ Rⁿ, kde n = 10, m = 1. Matice A obsahuje namˇeˇrené napˇet´ı, vektor pozorován´ı b namˇeˇrené proudy, a neznámá x je vodivost. Body jsou proloˇzeny lineárn´ı funkc´ı. Spojnice mezi body a funkc´ı ukazuj´ı odchylky mˇeˇren´ı, které se minimalizuj´ı (resp. souˇcet jejich kvadrát˚u) pˇri lineárn´ı regresi.

2.2 Zaveden´ı a odvozen´ı ˇ reˇ sen´ı pomoc´ı SVD

Nyn´ı se budeme zabývat úlohou kde obˇe sloˇzky namˇeˇrených dat obsahuj´ı chybu, tedy budeme opravovat oboje U i I, viz obrázek 2.2.

K ˇreˇsen´ı této úlohy potˇrebujeme následuj´ıc´ı definici.

Definice 4 ( Úplný problém nejmenˇs´ıch ˇctverc˚u). Necht’ existuj´ı A ∈ R^n×ma b ∈ Rⁿ. Upln´´ ym problémem nejmenˇs´ıch ˇctverc˚u nazveme úlohu nalezen´ı nejmenˇs´ı opravy matice A a vektoru pravé strany tak, aby opravená úloha byla kompatibiln´ı, t.j.

min k[r, E]k_F tak, ˇze (b + r) ∈ R(A + E).

Pak existuje vektor x ∈ R^m, který ˇreˇs´ı opravenou úlohu. Problém lze tedy také

(28)

0 2 4 6 8 10 12 14 16 18 0

2 4 6 8 10 12 14 16 18

U [V]

I [A]

Obrázek 2.2: Stejná data jako na obrázku 2.1 jsou nyn´ı proloˇzena lineárn´ı funkc´ı tak, ˇze se minimalizuje vzdálenost bod˚u mˇeˇrená kolmo k této funkci (resp. souˇcet kvadrát˚u tˇechto vzdálenost´ı). Tento pˇr´ıstup se nazývá ortogonáln´ı regrese.

formulovat n´asleduj´ıc´ım zp˚usobem

min k[r, E]k_F tak, ˇze (A + E)x = b + r, x ∈ R^m. (2.2) V praxi se tak´e m˚uˇzeme setkat s oznaˇcen´ım ortogon´aln´ı regrese.

Základn´ı uˇcebnice týkaj´ıc´ı se úplného problému nejmenˇs´ıch ˇctverc˚u je [18]. Pokud se vrát´ıme k aproximaˇcn´ı úloze Ax ≈ b (2.1), kde matice A je model, b je vektor pozorován´ı a x je vektor hledaných dat. Pak budeme opravovat matici modelu a vektor pozorován´ı. Pro jednoduchost pˇredpokládejme, ˇze matice A ∈ R^n×m, n > m, má plnou sloupcovou hodnost, tedy

rank(A) = m.

Dále pro vektor b ∈ Rⁿplat´ı z aproximaˇcn´ı úlohy podm´ınka b /∈ R(A), pak hodnost matice A rozˇs´ıˇrené o sloupec vektoru b vzroste o jedna, tedy

rank([b, A]) = m + 1.

(29)

Pomoc´ı SVD lze matici [b, A] zapsat v dyadick´em tvaru

[b, A] = SΘW^T =

m+1

X

j=1

θ_js_jw_j^T, (2.3)

pro kter´y plat´ı

θ₁ ≥ θ₂ ≥ . . . ≥ θ_m ≥ θ_m+1 > 0. (2.4) Nyn´ı budeme hledat nejmenˇs´ı opravu [r, E], tak aby platilo

rank([b, A] − [r, E]) < m + 1.

Podle vˇety 7 bude oprava tvaru

[r, E] = −θ_m+1s_m+1w^T_m+1. Pak opraven´a matice [b, A] bude tvaru

[bb, bA] =

m

X

j=1

θ_js_jw_j^T.

V dalˇs´ı sekci se budeme zabývat, zda opravená soustava bAx = bb má ˇreˇsen´ı.

2.3 Existence a jednoznaˇ cnost ˇ reˇ sen´ı

V této sekci se budeme zabývat zda opravená soustava bAx = bb má ˇreˇsen´ı a zda je toto ˇreˇsen´ı jednoznaˇcné, tedy zda má p˚uvodn´ı úloha Ax ≈ b ˇreˇsen´ı ve smyslu TLS.

Vrat’me se tedy k úloze Ax ≈ b. Jednoduchými úpravami lze úlohu zapsat ve tvaru [b, A] −1

x

≈ 0.

Obdobnˇe pro opravenou soustavu dostaneme [(b + r), (A + E)] −1

x

= 0 respektive

[bb, bA] −1 x

= 0. (2.5)

Zapiˇsme matici [bb, bA] v jej´ım dyadick´em tvaru

[bb, bA] =

m

X

j=1

θjsjw_j^T.

(30)

Nyn´ı obˇe strany vyn´asob´ıme vektorem w_m+1 a z´ısk´ame [bb, bA]w_m+1 =

m

X

j=1

θ_js_jw_j^Tw_m+1 =

m

X

j=1

θ_js_jhw_m+1, w_ji .

Protoˇze sloupce matice W , tedy vektory wj, jsou navz´ajem ortonorm´aln´ı, plat´ı [bb, bA]w_m+1 = 0.

Uvaˇzujme vektor w_m+1 ve tvaru

w_m+1 =





 ν₁ ν₂ ... ν_m+1





 .

Pokud ν₁ 6= 0 potom m˚uˇzeme zavést pomocný vektor w⁰_m+1, který bude tvaru

w_m+1⁰ = −1

ν₁ w_m+1 =







−1

−ν₂/ν₁ ...

−ν_m+1/ν₁





 .

Z porovnán´ı s (2.5) vyplývá, ˇze pokud ν₁ 6= 0, pak p˚uvodn´ı úloha Ax ≈ b má ˇreˇsen´ı ve smyslu TLS ve tvaru

x = −1 ν₁





 ν₂

... ν_m+1





=







−ν₂/ν₁ ...

−ν_m+1/ν₁





. Pokud ν₁ = 0 pak v´yˇse uveden´y postup zˇrejmˇe nelze pouˇz´ıt.

Obecnˇe se m˚uˇze stát, ˇze úloha má v´ıce neˇz jedno (resp. nekoneˇcnˇe mnoho) ˇreˇsen´ı.

Staˇc´ı si uvˇedomit, ˇze m˚uˇze nastat situace θ_m = θ_m+1 ve (2.4) (viz poznámka na konci kapitoly 1, str. 24) a tedy lze pro konstrukci ˇreˇsen´ı pouˇz´ıt jiný pravý singulárn´ı vektor. M˚uˇze se ale také stát, ˇze ˇreˇsen´ı v˚ubec nemá. Následuj´ıc´ı vˇeta formuluje postaˇcuj´ıc´ı podm´ınku pro existenci a jednoznaˇcnost ˇreˇsen´ı TLS.

Vˇeta 8 (Golub, Van Loan). Uvaˇzujme matici A ∈ R^n×m, pro kterou plat´ı n > m a rank(A) = m, t.j. A má lineárnˇe nezávislé sloupce. Dále necht’ plat´ı, ˇze b /∈ R(A).

Uvaˇzujme singul´arn´ı rozklady

A = U ΣV^T z (1.14)–(1.15) a

[b, A] = SΘW^T z (2.3)–(2.4).

Jestliˇze

σ_m > θ_m+1, (2.6)

pak pro matici [b, A] plat´ı, ˇze ν1 6= 0 a TLS problém má jednoznaˇcné ˇreˇsen´ı.

(31)

D˚ukaz viz [8] pˇr´ıpadnˇe [18]. Zde lze také nalézt diskusi týkaj´ıc´ı se pˇr´ıpadu s v´ıce neˇz jedn´ım ˇreˇsen´ım, pak hledáme ˇreˇsen´ı minimáln´ı v normˇe. Je vhodné zd˚uraznit, ˇze výˇse uvedená vˇeta formuluje pouze postaˇcuj´ıc´ı podm´ınku nikoli nutnou.

(32)

3 Core probl´ em

V pˇredchoz´ı kapitole jsme se zabývali ˇreˇsen´ım úlohy Ax ≈ b, respektive ˇreˇsen´ım opravené úlohy soustavy (A + E)x = b + r. Vyslovili jsme také postaˇcuj´ıc´ı podm´ınku pro existenci a vlastnˇe také jednoznaˇcnost ˇreˇsen´ı úlohy pomoci TLS. Zjistili jsme také, ˇze ne vˇzdy existuje ˇreˇsen´ı opravené úlohy pomoc´ı TLS. Nyn´ı se budeme zabývat metodou, která pro danou úlohu i jej´ı opravu vˇzdy nalezne ˇreˇsen´ı, tedy pokud úloha má ˇreˇsen´ı.

3.1 Ortogon´ aln´ı transformace ´ ulohy

Uvaˇzujme opˇet ´ulohu (2.1)

Ax ≈ b, kde A ∈ R^n×m, x ∈ R^m, b ∈ Rⁿ. (3.1) Dále zaved’me ortogonáln´ı matice P a Q odpov´ıdaj´ıc´ıch velikost´ı. Nyn´ı rovnici (3.1) vynásob´ıme zprava P^T a, s vyuˇzit´ım toho, ˇze Q je ortogonáln´ı matice, t.j. QQ^T = I, dostaneme

(P^TAQ)(Q^Tx) ≈ P^Tb. (3.2)

Zkr´acenˇe tak´e m˚uˇzeme zapsat

Aex ≈ ee b, kde A = Pe ^TAQ, ex = Q^Tx, eb = P^Tb.

Pouˇzijme stejnou ortogonáln´ı transformaci na opravenou úlohu (2.2), kde matice A a vektor b jsou opravené matic´ı E respektive vektorem r. Dostaneme

(P^T(A + E)Q)(Q^Tx) = (P^T(b + r)). (3.3) Protoˇze Frobeniova norma je ortogon´alnˇe invariantn´ı [2], plat´ı

min

[P^Tr, P^TEQ]

F = min

P^T[r, EQ]

F = min k[r, E]k_F. Rovnici (3.3) m˚uˇzeme tak´e zapsat ve tvaru

( eA + eE)ex = eb +r,e

kde eE = P^TEQ aer = P^Tr. Tedy hledan´y vektor x bude roven

x = Qex. (3.4)

(33)

Pˇredpokl´adejme, ˇze matice P a Q transformuj´ı rozˇs´ıˇrenou matici [eb, eA] do tvaru [eb, eA] = P^T[b, A] 1 0

0 Q

= b₁ A11 0 0 0 A₂₂

. (3.5)

Ulohu (3.2) pak m˚´ uˇzeme pˇrepsat ve tvaru

b₁ A₁₁ 0 0 0 A₂₂





−1 x₁ x₂



≈ 0 neboli, po pron´asoben´ı

−b₁+ A₁₁x₁ A₂₂x₂

≈ 0.

P˚uvodn´ı úloha se t´ım rozpadá na dva nezávislé podproblémy A₁₁x₁ ≈ b₁ a A₂₂x₂ ≈ 0.

Vid´ıme, ˇze druhý podproblém má ˇreˇsen´ı v klasickém slova smyslu (jinými slovy, je to soustava rovnic, nikoliv aproximaˇcn´ı úloha v pravém slova smyslu). Protoˇze zpravidla hledáme ˇreˇsen´ı minimáln´ı v normˇe (nemáme-li nˇejakou dodateˇcnou informaci k úloze a dobrý d˚uvod volit ˇreˇsen´ı jiné neˇz minimáln´ı), je nasnadˇe, ˇze jako ˇreˇsen´ı druhého podproblému zvol´ıme x₂ = 0. Druhý podproblém tedy nemá vliv na ˇreˇsen´ı p˚uvodn´ı úlohy. Zbývá vyˇreˇsit prvn´ı podproblém. Pokud x₁ bude ˇreˇsen´ım prvn´ıho podproblému, pak z (3.4) vid´ıme, ˇze

x = Qex = Q x₁ x₂

= Q x₁ 0

bude ˇreˇsen´ım p˚uvodn´ı ´ulohy.

3.2 Dovˇ etek k vˇ etˇ e 8

Oznaˇcme

Σ(M ) a min Σ(M )

mnoˇzinu vˇsech singulárn´ıch ˇc´ısel matice M vˇcetnˇe násobnost´ı, resp. nejmenˇs´ı sin- gulárn´ı ˇc´ıslo matice M , t.j.

|Σ(M )| = rank(M ).

Ze vztahu (3.5) zˇrejmˇe plat´ı

Σ(A) = Σ A₁₁ 0 0 A₂₂

= Σ(A11) ∪ Σ(A22), (3.6) proto, ˇze P a Q jsou ortogonáln´ı matice a ˇze mnoˇzina singulárn´ıch ˇc´ısel blokovˇe diagonáln´ı matice je sjednocen´ım mnoˇzin singulárn´ıch ˇc´ısel jejich diagonáln´ıch blok˚u.

Tak´e ale plat´ı

Σ([b, A]) = Σ b₁ A₁₁ 0 0 0 A₂₂

= Σ([b₁, A₁₁]) ∪ Σ(A₂₂). (3.7)

(34)

V pˇredpokladech vˇety 8 je, ˇze A má m´ıt lineárnˇe nezávislé sloupce. To je splnˇeno tehdy a jen tehdy kdyˇz obˇe matice A₁₁ a A₂₂ maj´ı lineárnˇe nezávislé sloupce. Dále b 6∈ R(A) tehdy a jen tehdy kdyˇz b₁ 6∈ R(A₁₁). Ze vztah˚u (3.6) a (3.7) plyne následuj´ıc´ı d˚uleˇzitý d˚usledek.

Necht’ A₁₁ a A₂₂ maj´ı pln´y sloupcov´y rank a b₁ 6∈ R(A₁₁). Pokud min Σ(A₂₂) < min{ min Σ(A₁₁) , min Σ([b₁, A₁₁)},

pak min Σ(A₂₂) = min Σ(A) = min Σ([b, A]). Vid´ıme, ˇze aˇckoliv jsou splnˇeny pˇredpoklady vˇety 8, nerovnost (2.6) nenastane. Nejmenˇs´ı singul´arn´ı ˇc´ısla

matic A a [b, A] jsou identick´a.

Porovnejme právˇe zformulovaný d˚usledek s tvrzen´ımi uvedenými v úvodu textu, viz [17], [5], resp. [6]. Poznamenejme jeˇstˇe, ˇze d´ıky tzv. prokládán´ı singulárn´ıch ˇc´ısel plat´ı min Σ(A₁₁) ≥ min Σ([b₁, A₁₁]), viz [14]. V pˇr´ıpadˇe core problému dokonce plat´ı min Σ()A₁₁ > min Σ([b₁, A₁₁]), viz [12], viz také vˇeta 9 v následuj´ıc´ı sekci.

3.3 Zaveden´ı core probl´ emu

V pˇredchoz´ı sekci jsme ukázali jak lze úlohy (3.1) rozloˇzit na dva podproblémy A₁₁x₁ ≈ b₁ a A₂₂x₂ ≈ 0.

V ˇclánku [12] bylo ukázáno, ˇze ortogonáln´ı transformace ve tvaru (3.2), vedouc´ı na bloky tvaru (3.5) vˇzdy existuje (m˚uˇze se ovˇsem stát, ˇze matice A₂₂má nulový poˇcet ˇrádk˚u nebo sloupc˚u). Nav´ıc vˇzdy existuje taková, ˇze matice A₁₁ má minimáln´ı a A₂₂ maximáln´ı dimenzi. Vzhledem k d˚uleˇzitosti tohoto pojmu zavedeme jeho formáln´ı definici.

Definice 5 (Core probl´em). ˇRekneme, ˇze A₁₁x₁ ≈ b₁ je core probl´em v aproximaˇcn´ı

´

uloze (3.1) za pˇredpokladu, ˇze [b₁, A₁₁] má minimáln´ı dimenzi (a matice A₂₂ ma- ximáln´ı dimenzi).

Nyn´ı uved’me vˇetu zab´yvaj´ıc´ı se vlastnostmi core probl´emu.

Vˇeta 9. Necht’ Ax ≈ b, b /∈ R(A) a A11x1 ≈ b1, kde A11 ∈ R^n×m, x1 ∈ R^m, b1 ∈ Rⁿ a b₁ ∈ R(A/ ₁₁) je core probl´em pak

• Matice A₁₁ m´a plnou sloupcovou hodnost, tedy rank(A₁₁) = m.

Protoˇze b₁ ∈ R(A/ ₁₁), pak n > m.

(35)

• Rozˇs´ıˇrená matice [b₁, A₁₁] ∈ R^n×(m+1) má plnou ˇrádkovou hodnost, tedy rank([b₁, A₁₁]) = n,

a tedy n = m + 1.

• Obˇe matice A11 a [b1, A11] maj´ı jednoduch´a singul´arn´ı ˇc´ısla σ₁⁰ > σ⁰₂ > . . . > σ⁰_m > 0,

θ⁰₁ > θ⁰₂ > . . . > θ⁰_m+1 > 0, pro kter´a plat´ı

θ₁⁰ > σ₁⁰ > θ⁰₂ > σ⁰₂ > . . . > σ⁰_m > θ_m+1⁰ .

Z právˇe uvedené vˇety z´ıskáváme d˚uleˇzitou vlastnost core problému nezávisle na p˚uvodn´ı úloze Ax ≈ b tedy, ˇze core problém vˇzdy splˇnuje pˇredpoklady vˇety 8.

Podle vˇety 8 má core problém A₁₁x₁ ≈ b₁ ˇreˇsen´ı a toto ˇreˇsen´ı je právˇe jedno nezávisle na p˚uvodn´ı úloze Ax ≈ b.

Je d˚uleˇzité upozornit, ˇze matice A₁₁ a vektor b₁ nejsou dány jednoznaˇcnˇe. Zˇrejmˇe pro libovolné ortogonáln´ı matice P₁ a Q₁ plat´ı, ˇze

(P₁^TA₁₁Q₁)(Q^T₁x₁) ≈ P₁^Tb₁ je opˇet core probl´em.

3.4 Tvary core probl´ emu

Core problém lze vyjádˇrit ve dvou základn´ıch tvarech tzv. SVD tvar a bidiagonáln´ı tvar.

3.4.1 Core probl´ em v SVD tvaru

V pˇredminulé sekci o ortogonáln´ı transformaci jsme pomoc´ı ortogonáln´ıch matic P a Q transformovali p˚uvodn´ı úlohu (2.1) do tvaru (3.2). Transformac´ı se úloha rozpadla na dva podproblémy a sice na podproblém A₁₁x₁ ≈ b₁ minimáln´ı dimenze nazvaný core problém a podproblém A₂₂x₂ ≈ 0 maximáln´ı dimenze. Uvaˇzujme obecný core problém A11x1 ≈ b1 a SVD matice A11

A₁₁ = U₁₁Σ₁₁V₁₁^T.

Vynásob´ıme-li rozˇs´ıˇrenou matici [b₁, A₁₁] zleva ortogonáln´ı matic´ı U₁₁ a zprava or- togonáln´ı matic´ı diag(1, V₁₁)^T, dostane tvar

U₁₁[b₁, A₁₁] 1 0 0 V₁₁

= [U₁₁^Tb₁|Σ₁₁].

V ˇclánku [12] bylo ukázáno, ˇze U₁₁^Tb₁ má vˇsechny sloˇzky nenulové. Tento zápis se nazývá SVD tvar core problému.

(36)

Definice 6 (SVD tvar core problému). Core problém se nazývá core problémem v SVD tvaru pokud

[b₁, A₁₁] =







β₁ σ⁰₁ ... . ..

σ_m⁰ β_m+1 0 . . . 0





 ,

kde σ₁⁰ > σ₂⁰ > . . . > σ_m⁰ > 0 a β_j 6= 0, j = 1, . . . , m + 1.

3.4.2 Bidiagon´ aln´ı tvar core probl´ emu

Nyn´ı se budeme zabývat bidiagonáln´ım tvarem core problému. Pomoc´ı ortogonáln´ıch matic P , Q lze transformovat matici [b, A] do matice bidiagonáln´ıho tvaru pomoc´ı napˇr´ıklad Householderovy reflexe. Transformace pomoc´ı Householderových reflex´ı spoˇc´ıvá ve vytvoˇren´ı posloupnosti reflex´ı H₁, H₂, . . . , které stˇr´ıdavˇe násob´ıme s matic´ı [b, A] zleva a zprava. Kde matice H_2j−1 nuluj´ı poddiagonáln´ı prvky v j-tém sloupci rozˇs´ıˇrené matice a matice H2j nuluj´ı naddiagonáln´ı prvky v j-tém ˇrádku systémové matice. Householderova reflexe bude vypadat

H_2j−1. . . H₁[b, A] eH₂. . . eH_2j kde matice eH_2j jsou tvaru

He_2j = 1 0 0 H_2j

.

Bidiagonalizaci ukonˇc´ıme ve chv´ıli kdy dojde v transformaci k oddˇelen´ı podproblém˚u, pak tedy matice A11bude v bidiagonáln´ım tvaru. Struktura matice A22nen´ı d˚uleˇzitá.

V ˇclánku [12] bylo ukázáno, ˇze podproblém oddˇelený bidiagonalizac´ı je opˇet core problémem, tentokrát v bidiagonáln´ım tvaru.

Definice 7 (Bidiagonáln´ı tvar core problému). Core problém se nazývá core pro- blémem v bidiagonáln´ım tvaru pokud

[b₁, A₁₁] =







δ₁ γ₁ δ₂ γ₂

. .. ...

δm γm

δ_m+1





 ,

kde γ` > 0, ` = 1, . . . , m, a δ_j > 0, j = 1, . . . , m + 1.

(37)

4 Automatick´ e generov´ an´ı aproximaˇ cn´ıch

´

uloh s core probl´ emem s pˇ redem zn´ am´ ymi parametry

V této kapitole si pop´ıˇseme software, který generuje core problémy A₁₁x₁ ≈ b₁ v SVD tvaru i v bidiagonáln´ım tvaru podle pˇredem zadaných parametr˚u, resp. aproximaˇcn´ı

´

ulohy Ax ≈ b s core problémem s pˇredem známou velikost´ı a s pˇredem známými parametry. Software jsme vytvoˇrili v prostˇred´ı programu Matlab.

4.1 Program cpSVD. Core probl´ em v SVD tvaru

C´ılem tohoto programu bude vytvoˇrit core probl´em, tedy matici A₁₁a pravou stranu b₁, v SVD tvaru. Tedy

[b₁, A₁₁] =







β₁ σ⁰₁ ... . ..

σ_m⁰ β_m+1 0 . . . 0





 ,

kde σ₁⁰ > σ₂⁰ > . . . > σ_m⁰ > 0 a β_j 6= 0, j = 1, . . . , m + 1.

Vstupn´ımi parametry programu cpSVD budou obecnˇe velikost vytvoˇreného core problému aproximaˇcn´ı úlohy a typ vytvoˇren´ı SVD tvaru tedy zp˚usob urˇcen´ı hodnot βj a σj. U r˚uzných typ˚u pak m˚uˇzou pˇrib´ıt dalˇs´ı parametry, které ovlivn´ı urˇcen´ı hodnot β_j a σ_j.

Program cpSVD. Zad´an´ı rozmˇeru a typ core probl´emu function [b_1,A_11] = cpSVD(m,typ,par)

% m - velikost core probl´emu

% typ = 11 - náhodná ˇc´ısla beta a sigma s rovnomˇerným

% rozdˇelen´ım z intervalu (0,1); pˇr´ıkaz rand

% typ = 12 - náhodná ˇc´ısla beta a sigma s normáln´ım

% rozdˇelen´ım N(0,1); pˇr´ıkaz randn