Современные проблемы прикладной математики. Часть 2. Практикум

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ

ТОМСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ СИСТЕМ УПРАВЛЕНИЯ И РАДИОЭЛЕКТРОНИКИ (ТУСУР)

Ю.Е. Воскобойников А.А. Мицель

СОВРЕМЕННЫЕ ПРОБЛЕМЫ ПРИКЛАДНОЙ МАТЕМАТИКИ

Часть 2. Практикум

Учебное пособие

ТОМСК 2016

УДК 519.2 ББК 22.172

В650

Воскобойников Ю. Е., Мицель А.А.

Современные проблемы прикладной математики. Часть 2. Практикум: учебное пособие/ Ю. Е. Воскобойников, А.А. Мицель/ Томский гос. ун-т систем управления и радиоэлектроники (ТУСУР). – Томск, 2016. – 52с.

В учебном пособии в первой части приводится системное изложение одного из разделов прикладной математики, связанного с устойчивыми методами и алгоритмами решения систем линейных алгебраических уравнений, возникающих при параметрической идентификации моделей. Основное внимание уделяется построению решений с минимальной ошибкой или с требуемыми точностными характеристиками, а также учету имеющейся априорной информации об искомом решении. Во второй части приводится описание практических занятий по созданию алгоритмов построения нормального псевдорешения и регуляризированных решений систем линейных алгебраических уравнений.

Учебное пособие предназначено для магистрантов направления «Прикладная математика и информатика». Результаты будут полезны также широкому кругу студентов, магистрантов, аспирантов, исследователей, занимающихся решением задач параметрической идентификации и обработки экспериментальных данных.

1

2

ОГЛАВЛЕНИЕ
ПРАКТИЧЕСКОЕ ЗАНЯТИЕ №1
Построение нормального псевдорешения СЛАУ .	5
§ 1.1. Постановка задачи ………………..……………..…	5
§ 1.2. SVD-алгоритм построения нормального
псевдорешения…………..….…………………………….	6
Задание 1.1…………………………………………………	8

Задание 1.2…………………………………………………. 12

ПРАКТИЧЕСКОЕ ЗАНЯТИЕ №2

Построение регуляризованного решения СЛАУ . 14

§2.1. Байесовский регуляризирующий алгоритм………... 14 Задание 2.1…………………………………………………. 16 §2.2. Оптимальный регуляризирующий SVD-алгоритм 17 Задание 2.2…………………………………………………. 20 §2.3. Построение регуляризованного решения при неполной информации…………………………………… 21 Задание 2.3………………………………………………… 26 §2.4. Регуляризирующий SVD-алгоритм………………… 27 Задание 2.4…………………………………………………. 30

ПРАКТИЧЕСКОЕ ЗАНЯТИЕ №3

Алгоритмы выбора параметра регуляризации . 35

§3.1. Выбор параметра регуляризации на основе критерия оптимальности…………………………………. 35

Задание 3.1………………………………………………… 36 §3.2. Алгоритм выбора параметра регуляризации с использованием SVD-разложения на основе критерия оптимальности…………………………………………… 36 Задание 3.2……………………………………………….. 38 §3.3. Алгоритм выбора параметра регуляризации основе статистического принципа невязки……………………. 38

Задание 3.3………………………………………………... 40

§3.4. Алгоритм поиска αV с использованием SVD разложения……………………………………………….. 41 Задание 3.4………………………………………………. 42

ПРАКТИЧЕСКОЕ ЗАНЯТИЕ №4

Локальная регуляризация ………………………….. 43

§4.1. Векторный параметр регуляризации……………. 43 Задание 4.1………………………………………………… 45

ПРАКТИЧЕСКОЕ ЗАНЯТИЕ №5

Дескриптивные алгоритмы решения СЛАУ …… 46

§5.1. Глобальный дескриптивный регуляризирующий алгоритм ……………. …………….……………………… 46 Задание 5.1………………………………………………… 49 § 5.2. Локальный дескриптивный регуляризирующий алгоритм ……………. …………….………………… 49 Задание 5.2…………………………………………………. 51

ЛИТЕРАТУРА……………………………………………….… 52

3

4

Практическое занятие №1

Построение нормального псевдорешения СЛАУ

§1.1. Постановка задачи

Дана система линейных алгебраических уравнений

матричном виде
					Kϕ = f ,	(1.1)
где	K				– матрица размером N × M (N строк и M столбцов),
ϕ	–		вектор размерности M (содержит M проекций),			f –
вектор размерности N
fɶ =					+η
		f
Здесь					- вектор точной правой части, η - вектор ошибок.
			f
			Предположим, что матрица K имеет размеры N × M .

Вектор ϕHK размерностью M называют псевдорешением (или

решением МНК), если он доставляет минимум следующему функционалу

ΨHK (ϕ) =				f − Kϕ				2 = ( f − Kϕ)T ( f − Kϕ)	(1.2)

среди всех векторов евклидова пространства EM .

Решение, обеспечивающее минимум функционалу (2), является решением следующей СЛАУ

KT Kϕ	HK	= KT f ,	(1.3)
которая называется системой нормальных уравнений.
В отличие от исходной системы Kϕ = f			эта система
всегда разрешима, т.е. для любой правой части f			существует

псевдорешение ϕHK . Если матрица	K имеет ранг, равный M ,
то
ϕHK = (KT K)−1 KT f .		(1.4)
Сингулярным разложением	прямоугольной	N × M
матрицы K (коротко: SVD-разложением) называется
представление:
K = UΛV T ,		(1.5)

где U – ортогональная ( N × N )-матрица, V – ортогональная ( M × M )-матрица, Λ – ( N × M )-матрица вида

λ 1 0 0 0

0λ 2 0 0

	Λ	=	0	0	λ 3 0			,	(1.6)
			0	0	0	λ M
			0	0	0	0
в которой последние N – M строки содержат только нулевые
элементы.	Величины				λj	≥ 0,	j = 1,...,M ,		называются
сингулярными числами матрицы K, и в дальнейшим полагаем, что
λj упорядочены по убыванию, т.е. λj								≥ λj+1 . Напомним, что

матрица B называется ортогональной, если имеет место тождество

BT B = BBT = I

§1.2. SVD-алгоритм построения нормального
псевдорешения
Введем векторы
y = UT f , x =V T ϕ	(1.7)

5

6

размерностью N и M соответственно. Тогда с учетом (1.5)

систему Kϕ = f можно преобразовать к эквивалентной системе:

	λj xj = yj ,	j = 1,...,M;		(1.8)
	0 = yj ,	j = M +1,...,	N,

которая хорошо характеризует «информативность» правой части:

чем меньше сингулярное число λj	, тем с меньшим весом проекция
xj входит в правую часть.	Предельный				случай λj = 0,
p +1≤ j ≤ M , говорит о вырожденности K.
		N
Очевидно, что невыполнение условия		∑		yj		= 0 говорит о
несовместности исходной системы.		j=M +1

Сучетом ортогональности матриц U, V и соотношений

(5)функционал (2) можно записать в виде

	p
	ΨНК (ϕ) = ΨНК (x) = ∑( yi − λj xj )2 +
	j=1		.
	M	N
	+ ∑ (yj − 0 xj )2 +	∑ y2j
	j= p+1	j=M +1

Третье слагаемое обусловлено несовместностью исходной системы, и не зависит от x . Второе слагаемое отражает вырожденность системы, и, следуя определению нормального псевдорешения, проекции xj , входящие во второе слагаемое,

следует принять равным 0. Тогда минимум функционала достигается на векторе x+ размерности p с элементами

	x+j =	yj	, j = 1,..., p ,	(1.9)
		λj
	а нормальное псевдорешение ϕ+		выражается как

p
ϕ+ = ∑x+j vj .	(1.9)

j=1

Напомним, что λj > 0, j = 1,..., p , где p – ранг матрицы K .

Задание 1.1. Построить нормальное псевдорешение с помощью пакета Mathcad

Рассмотрим две функции Mathcad, которые потребуются

для	построения	нормального				псевдорешения	ɶ+
							ϕ	,
определяемого выражением
		ɶ+	p		ɶ
				П	f ,uj
		ϕ = ∑				vj ,	(1.10)
					λj
			j=1

где практический ранг pП матрицы K определяется количеством сингулярных чисел λj , удовлетворяющих условию:

λj	≥ γ 0 ,	(1.11)
λmax

где γ 0 – достаточно малая величина (10−10 10−8 ).

Функция svds. Обращение имеет вид svds(K). Вычисляет

вектор размерности M , состоящий из сингулярных чисел λj матрицы K, которые расположены в убывающем порядке.

Функция svd. Обращение имеет вид		svd(K).	Вычисляет
матрицу UV размером (N + M )× M . Первые			N строк
этой матрицы соответствуют матрице U		размером N × M ,

которая определяет первые M столбцов матрицы U , т.е.

7

8

U =		U		uM +1 uN		.	(1.12)

Последние M строк матрицы UV содержат матрицу V размером

M × M .

Заметим, что отсутствие в матрице U последних N − M

столбцов матрицы U обусловлено тем, что эти столбцы не участвуют в вычислении нормального псевдорешения и поэтому во многих программных реализациях SVD-разложения эти столбцы не вычисляются.

Функция submatrix. Обращение имеет вид submatrix (K, i1, i2, j1, j2). Формирует новую матрицу из элементов матрицы

K, стоящих с	i1 по i2 строках и с j1 по j2 столбцах матрицы K.
Пример	1. Дана матрица K размером 6×3. Необходимо

вычислить сингулярные числа и матрицы U,V .

Решение. На рис. 1.1 показан фрагмент документа Mathcad, выполняющий требуемые вычисления. Здесь же приведены вычисление числа обусловленности по формуле

cond(K) = λmax /λmin

и проверка ортогональности столбцов матриц U,V .☻

Перейдем к подпрограмме-функции (П-Ф) Ps_Solve, осуществляющей построение нормального псевдорешения СЛАУ по формуле (1.1). Обращение к П-Ф имеет вид:

Ps_Solve(K,f,γ 0 ).	(1.13)
Формальные параметры: K – матрица системы размером
N × M , f – правая часть системы, γ 0	– переменная
вещественного типа, входящая в условие (1.11).

На рис. 1.2 приведен фрагмент документа Mathcad с текстом П-Ф Ps_Solve.

Замечание 1. При обработке матриц в пакете Mathcad часто используется операция формирования вектора из определенного столбца матрицы. Для этого надо ввести имя матрицы, затем нажать клавиши [Ctrl+6] и в появившихся вверху угловых скобках задать нужный номер столбца. Например, в П-Ф

Ps_Solve стоят операции U j,V j . ♦

Рис. 1.1. Сингулярное разложение матрицы K

9

10

fɶ −

ɶ

f

−ϕ

ϕ

≤ cond(K)

.

ϕ

f

Действительно

1.103 103 ≤ cond(K) 3.232 10−3 = 4.609 103 .

Задание 1.2.

Вычислить нормальное решение с помощью обратной матрицы по формуле (1.4) для точной и зашумленной правой части и сравнить с решением, полученным SVD – алгоритмом.

Решение. Два обращения к Ps_Solve показаны на рис. 1.3.

Рис. 1.2. Текст подпрограммы-функции Ps_Solve

Пример 2. Матрица K размером 5×3 формируется с использованием П-Ф Form_K (фрагмент документа показан на рис. 1.3). Число обусловленности 1.426 106 . Для заданного вектора ϕ вычислены два вектора: вектор «точной» правой части f и вектор «зашумленной» правой части fη с

относительной погрешностью 3.232 10−3 (см. рис. 1.3). По этим двум векторам необходимо построить нормальные псевдорешения с использованием П-Ф Ps_Solve.

Здесь вычислены относительные ошибки двух псевдорешений: псевдорешение ϕ1ps , построенное по точной

правой части, и псевдорешение ϕ2ps , построенное по

искаженной правой части fη . Несмотря на маленькую погрешность исходных данных, относительная ошибка решенияϕ2ps достигает большой величины 1.102 103 , и эта ошибка удовлетворяет неравенству

11

12

Рис 1.3. Построение нормальных псевдорешений

Практическое занятие №2

Примечание. При наличии научной задачи магистранта,

Построение регуляризованного решения СЛАУ

которую можно свести к системе линейных уравнений,

построить нормальное псевдорешение для этой задачи.

§2.1. Байесовский регуляризирующий алгоритм

Исходной информацией являются:

•

априорное

распределение

p(ϕ )

искомого

вектора

решения

ϕ

является

нормальным,

т.е.

−

1

ϕ − mϕ

2

mϕ и

Vϕ -

p(ϕ) = const exp

2

−1

,

где

Vϕ

математическое ожидание и матрица ковариаций;

•

условное распределение p( fɶ

ϕ ) ,

характеризующее

распределение

вектора

измерений

fɶ

при фиксированном

векторе

ϕ

является

нормальным,

ɶ

−

1

ɶ

2

p( f

ϕ ) = const exp

2

f − Kϕ

−1

.

Vη

2 −1

Здесь

запись

z − mz

означает:

Vz

z − mz V2z−1= (z − mz )T Vz−1 (z − mz ) .

Байесовское регуляризированное решение находится из системы линейных алгебраических уравнений

(KTV

−1K +V −1 )ϕ

Б

= KTV −1 fɶ +V −1m .

(2.1)

η

ϕ

η

ϕ

ϕ

Заметим, что вектор ϕБ

для указанных

распределений и

квадратичной

функции

потерь

П (ϕT ,ϕ ) =

ϕT −ϕ

2

максимизирует

значения

апостериорной

плотности

распределения

p(ϕ | fɶ) , которая имеет вид

13

14

ɶ

1

ɶ

2

1

2

p(ϕ | f ) = const exp

−

2

f

− Kϕ

−1

−

2

ϕ − mϕ

−1

,

Vη

Vϕ

а, следовательно, доставляет минимум функционалу

F	[ϕ] =		fɶ − Kϕ		2	−1	+		2 −1	.	(2.2)
								ϕ − m
Б					V			ϕ	V
					η				ϕ
Рассмотрим	точность						построенного				байесовского

регуляризированного решения при следующих предположениях:

M [η] = 0 ;	M ϕηT		= 0 .	(2.3)
Последнее условие означает, что проекция ϕ j				и ηi не

коррелированны между собой. Определим вектор ошибки решения как

εБ = ϕБ −ϕ		.	(2.4)

Ошибка может быть представлена в следующем виде:

εБ = (KTVη−1K +Vϕ−1 )−1 (KTVη−1 (Kϕ +η) +Vϕ−1mϕ ) −

		= (KTV −1K +V −1 )−1		KTV −1η +V −1			(m			)	(2.5)
−ϕ								−ϕ			.
		η	ϕ		η	ϕ	ϕ

Математическое ожидание M [εБ ] равно:

M [εБ ] = bБ = (KTVη−1K +Vϕ−1 )−1Vϕ−1 (mϕ − Mϕ [ϕ]) . (2.6)

Рассмотрим случайную ошибку	ξБ = εБ − М [εБ ] ,	которая
представима в виде:
ξБ = (KTVη−1K +Vϕ−1 )−1 KTVη−1η ,		(2.7)
с математическим ожиданием	M [ξБ ] = 0 и	матрицей
ковариации

Vξ	= (KTVη−1K +Vϕ−1 )−1	KTVη−1K (KTVη−1K +Vϕ−1 )−1 .	(2.8)
	Б

Вектор «полной» ошибки байесовского решения можно записать в виде:

∆	Б j	= b	+ V , j = 1,...,M .			(2.9)
		Б j		ξ jj

Задание 2.1.

1.Задать матрицу системы K , вектор решения ϕ , среднее значение mϕ и матрицу ковариации Vϕ .

2.Для заданной матрицы K и вектора решения ϕ вычислить

правую часть f .

3. Внести шум в правую часть, т.е. вычислить вектор

f = f + η, где η - вектор ошибки, вычисляемый с помощью нормального датчика случайных чисел для двух случаев:

а) ошибки не коррелированны. Значения СКО задать равными σηi = δi fi , i =1,..., N ; δi = 0.01,0.05,0.01.

Матрица ковариации в этом случае будет диагональной

Vη = diag (σ2η1 ...,σ2ηN ).

б) ошибки коррелированны. В этом случае вектор ошибок вычисляется по формуле

η =Vη1/2ξ,

где Vη - заданная матрица ковариации; ξ - вектор, проекции

которого состоят из нормальных случайных величин с единичной дисперсией и нулевым средним.

4. Вычислить решение по формуле (2.1) и вектор ошибки (2.9).

15

16

Примечание. При наличии научной задачи магистранта, которую можно свести к системе линейных уравнений, провести исследования этой задачи.

§2.2. Оптимальный регуляризирующий SVD-алгоритм

Регуляризированное решение запишем в виде

ϕ

R

=VRUTC−12

fɶ

(2.10)

η

или в скалярном виде

p

2

ϕR = ∑

λj

j+ +

λj

uj

,n vj ,

x

(2.11)

2

2

j=1

λj + qj

λj

+ qj

где матрица R размером M × N имеет следующую структуру:

r1

0

0

0

0

r2

0

0

R =

0

,

0

rM

0

r

=

λj

j =1,2,...,M

j

λ2

+ q

,

.

(2.12)

j

j

Величины rj

назовем регуляризирующими множителями.

Здесь n = C

−12η ,

+ =

v

+

+

x

,ϕ

,

- нормальное

точное

j

ϕ

η

j

псевдорешение,

построенное

для точной

правой части и

			pП	f ,uj
			+ = ∑		vj ; vj ,uj – j -й
	определяемое формулой ϕ
			j=1	λj

столбец матрицы V и U соответственно; p – ранг матрицы K.

Матрица ковариации Vη представлена в виде

V = σ 2	C	(2.13)
η η	η

Очевидно, от величины регуляризирующих множителей зависит ошибка решения ϕR , которую определим функционалом

∆(R) = Mη

+

2

ϕR −ϕ

,

(2.14)

которая может быть записана в виде

p

qj

2

p

λj

2

∆(R) = ∑

(

j+ )2

+ση2

∑

, (2.15)

x

λ

2

λ

2

j

+ q

j

+ q

j=1

j

j=1

j

Определим элементы матрицы R построенного

алгоритма путем нахождения

величин

qj ,

j =1,2,..., p, из

условия минимума функционала ∆(R) . Дифференцируя (2.15) по qj и приравнивая производную нулю, получаем оптимальное значение

=		σ 2
			η
qоптj				, j =1,2,..., p .	(2.16)
	(		j+ )2
		x
Оптимальное решение	ϕопт , построение при qj = qопт j				имеет

следующее SVD-представление:

17

18

p

λj

−1

ϕопт

= ∑

2 fɶ vj ,

uj

,Cη

(2.17)

2

+ q

j=1

λ

j

оптj

или в матричном виде:

−1

2 fɶ ,

ϕопт = Vp Rpопт U pT Cη

(2.18)

где Vp – матрица размера

M × p , составленная из p

первых

столбцов матрицы

V ;

Up

–

матрица

размера

N × p ,

составленная из

p

первых столбцов

матрицы

U ; Rpопт –

диагональная матрица размера p × p следующей структуры:

λ1

0

0

λ2

+ q

1

опт

1

0

λ1

0

R

=

λ22 + qопт

. (2.19)

pопт

2

0

0

0

λp

λ2 + q

p

оптp

Точность построенного решения определим вектором ошибки

εопт = ϕопт −ϕ +

оптимального решения ϕопт	в виде следующей суммы
векторов:

εопт = ϕопт −ϕопт +ϕопт −ϕ + = ξопт + bопт . (2.20)

Вектор

b

+

определяет систематическую

ошибку

= ϕ

опт

−ϕ

опт

решения ϕопт . При этом

bопт = Mη [εопт ] и поэтому вектор

b

можно

назвать

смещением

решения ϕопт .

Вектор

опт

опт

является случайным вектором с нулевым

ξопт = ϕопт −ϕ

средним

M [ξопт ] = 0

и характеризует случайную

ошибку

решения ϕопт . Таким образом, для вектора смещения имеем

+ = VRUT C−1/ 2

−

+

= V (RΛ − I)

+ .(2.21)

b

f

ϕ

x

= ϕ

опт

−ϕ

опт

Случайный вектор ξопт

ξопт

опт

= ϕопт = Vp Rpопт U pT Cη−12η

(2.22)

= ϕопт − ϕ

имеет нулевое среднее и матрицу ковариации

V = σ2V

R2

V

T .

(2.23)

ξ

η

p

pопт

p

Таким образом, полная ошибка равна

∆оптj

= bоптj

+

j = 1,..., p .

Vξ jj ,

(2.24)

Задание 2.2.

1.Задать матрицу системы K , вектор решения ϕ , среднее значение mϕ и матрицу ковариации Vϕ .

2.Для заданной матрицы K и вектора решения ϕ вычислить

правую часть f .

3. Внести шум в правую часть, т.е. вычислить вектор

f = f + η,

19

20