Сингулярное разложение матриц

@ZAK53 · Регистрация: 11.11.2013

Студворк — интернет-сервис помощи студентам

Не уверен, что это кому-нить нужно, но что делаю, о том пою.

Сингулярное разложение (SVD) - очень полезный, а иногда и необходимый аппарат для считальщиков. Применить можно много для чего, в частности - для решения линейных систем с плохо обусловленными или вообще неполного ранга матрицами (тестовый пример как раз с такой матрицей). Еще - для решения переопределенных (уравнений больше, чем неизвестных) систем методом Наименьших Квадратов.

В наш век, когда космические корабли бороздят и вычислительные мощности настольных ЭВМ огромны, вообще нет резона применять для решения систем более простые и экономные алгоритмы типа Гаусса. С ними можно получить "изумительные" результаты при плохой обусловленности.

Ниже привожу перевод на РВ ФОРТРАН 4 -подпрограммы из книги [Форсайт, Малькольм, Моулер. Машинные методы математических вычислений. М., 1982(?)], которая, в свою очередь, являлась переводом с небольшой модификацией АЛГОЛ-процедуры [Уилкинсон, Райнш. Справочник алгоритмов на АЛГОЛ].

Тестовый пример тоже взят из [Форсайт, Малькольм, Моулер].

Описание оставил фортрановское ... Нулевые индексы в массивах не используются, индексация начинается с 1. Потамушто. Единственное, что порадовало в РВ, так это то ли отсутствие оптимизатора, то ли его ненавязчивость. Благодаря этому проходят тесты на пренебрежимую малость типа
if A+anorm = anorm
и оптимизатор не превращает их в if A=0.

Черррт, исходник не проходит по размеру, попробую прицепить прицеп. Прицеп тоже не лезет, если кому интересно - возьмите с мого яндекс-диска: https://yadi.sk/d/wLnV4TYncg5iJ

@ZAK53 · 13.11.2014, 02:46 **[ТС]**

PureBasic
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
EnableExplicit
;---------------------------------------------------------------------------
 
Procedure.d maxD (a.d, b.d)
  Protected.d c
  If a => b : c=a : Else : c=b : EndIf
  ProcedureReturn c
EndProcedure
;-------------------------------------------------------------------------------------------
;    real*8 function SVD1(nm,m,n,w,matu,u,matv,v)
;
; подпрограмма вычисления сингулярного разложения матрицы A = UWV'
; действительной прямоугольной матрицы A(m,n),
; при этом используются двухдиагонализация посредством
; хаусхолдеровых отражений и вариант QR-алгоритма
;
;    integer     nm  >= max(m,n) !! На РВ этот параметр не нужен
;   * ,          m        ; число строк матрицы U
;   * ,          n        ; число столбцов матрицы U и порядок V
;    real*8      w(n)     ; n неотрицательных сингулярных чисел А
;    integer     matu     ; признак вычисления матрицы U
;    real*8      u(nm,n)  ; на входе содержит матрицу А, на выходе - U, т.е. исходная матрица не сохраняется.
;    integer     matv     ; признак вычисления матрицы V
;    real*8      v(nm,n)  ; ортогональная матрица
 
; примечание :    matu=1 , если нужно вычислять матрицу U
;                 из разложения , и matu=0 в противном
;                 случае
;                 matv=.1 , если нужно вычислять матрицу V
;                 из разложения , и matv=0 в противном
;                 случае
;
;  Возвращает 1 в случае успеха, и 0, если нет сходимости к сингулярному числу.
;
Procedure.i Svd1(m.i ,n.i ,Array W.d(1), matu.i, Array U.d(2), matv.i, Array V.d(2))
 
  Protected.i i,j,k,l,i1,k1,l1,mn,its
  Protected.d c,f,g,h,s,x,y,z,scale,anorm,eps
  Protected Dim Rv1.d(n)
;
; хаусхолдерово приведение к двухдиагональной форме
;
      g=0.0
      scale=0.0
      anorm=0.0
      For i=1 To n
         l=i+1
         rv1(i)=scale * g
         g=0.0
         s=0.0
         scale=0.0
         If i <= m
            For k=i To m
               scale=scale+Abs(u(k,i))
            Next
            If scale <> 0.0
               For k=i To m
                  u(k,i)=u(k,i)/scale     ; нормирование на SCALE
                  s=s+u(k,i) * u(k,i)
               Next
               f=u(i,i)
               g = -Sqr(s) * Sign(f)  ;g=-Sign(Sqrt(s),f)
               h=f * g-s
               u(i,i)=f-g
               If i <> n
                  For j=l To n
                     s=0.0
                     For k=i To m
                        s=s+u(k,i) * u(k,j)
                     Next
                     f=s/h
                     For k=i To m
                        u(k,j)=u(k,j)+f * u(k,i)
                     Next
                  Next
               EndIf
               For k=i To m
                  u(k,i)=scale * u(k,i)
               Next
            EndIf
         EndIf
         w(i)=scale * g
         g=0.0
         s=0.0
         scale=0.0
         If i <= m And i <> n
            For k=l To n
               scale=scale+Abs(u(i,k))
            Next
            If scale <> 0.0
               For k=l To n
                  u(i,k)=u(i,k)/scale
                  s=s+u(i,k) * u(i,k)
               Next
               f=u(i,l)
               g = -Sqr(s) * Sign(f)   ;g=-Sign(SQRT(s),f)
               h=f * g-s
               u(i,l)=f-g
               For k=l To n
                  rv1(k)=u(i,k)/h
               Next
               If i <> m
                  For j=l To m
                     s=0.0
                     For k=l To n
                        s=s+u(j,k) * u(i,k)
                     Next
                     For k=l To n
                        u(j,k)=u(j,k)+s * rv1(k)
                     Next
                  Next
               EndIf
               For k=l To n
                  u(i,k)=scale * u(i,k)
               Next
            EndIf
         EndIf
         anorm=MaxD(anorm,Abs(w(i))+Abs(rv1(i)))
      Next
;
; накопление правосторонних преобразований
;
      If matv
         For i=n To 1 Step -1
            If i <> n
               If g <> 0.0
                  For j=l To n
                     v(j,i)=(u(i,j)/u(i,l))/g
                  Next
                  For j=l To n
                     s=0.0
                     For k=l To n
                        s=s+u(i,k) * v(k,j)
                     Next
                     For k=l To n
                        v(k,j)=v(k,j)+s * v(k,i)
                     Next
                  Next
               EndIf
               For j=l To n
                  v(i,j)=0.0
                  v(j,i)=0.0
               Next
            EndIf
            v(i,i)=1.0
            g=rv1(i)
            l=i
         Next
      EndIf
;
; накопление левосторонних преобразований
;
      If matu
         mn=n
         If m < n :   mn=m : EndIf
         For i=mn To 1 Step -1
            l=i+1
            g=w(i)
            If i <> n
               For j=l To n
                  u(i,j)=0.0
               Next
            EndIf
            If g <> 0.0
               If i <> mn
                  For j=l To n
                     s=0.0
                     For k=l To m
                        s=s+u(k,i) * u(k,j)
                     Next
                     f=(s/u(i,i))/g
                     For k=i To m
                        u(k,j)=u(k,j)+f * u(k,i)
                     Next
                  Next
               EndIf
               For j=i To m
                  u(j,i)=u(j,i)/g
               Next
            Else
               For j=i To m
                  u(j,i)=0.0
               Next
            EndIf
            u(i,i)=u(i,i)+1.0
         Next
      EndIf
;
; диагонализация двухдиагональной формы
;
      For k=n To 1 Step -1
         k1=k-1
         its=0
         Repeat
           ; проверка возможности расщепления
           Repeat
             For l=k To 1 Step -1
                 l1=l-1
                 If Abs(rv1(l))+anorm = anorm : Break 2: EndIf
                 If (Abs(w(l1))+anorm = anorm) :  Break : EndIf
             Next
 
             c=0.0
             s=1.0
             For i=l To k
                f=s * rv1(i)
                rv1(i)=c * rv1(i)
                If Abs(f)+anorm = anorm  : Break 2: EndIf
                g=w(i)
                h=Sqr(f * f+g * g)
                w(i)=h
                c=g/h
                s=-f/h
                If Not matu : Continue : EndIf
                For j=1 To m
                   y=u(j,l1)
                   z=u(j,i)
                   u(j,l1)=y * c+z * s
                   u(j,i)=-y * s+z * c
                Next
             Next
             Break
           ForEver
 
; проверка сходимости
 
            z=w(k)
            If l = k  :  Break  : EndIf   ; выход из цикла Repeat
            If its = 100 : ProcedureReturn 0 : EndIf          ;нет сходимости к сингулярному числу
            its=its+1
            x=w(l)
            y=w(k1)
            g=rv1(k1)
            h=rv1(k)
            f=((y-z) * (y+z)+(g-h) * (g+h))/(2.0 * h * y)
            g=Sqr(f * f+1.0)
            f=((x-z) * (x+z)+h * (y/(f + g * Sign(f))-h))/x    ;   f=((x-z) * (x+z)+h * (y/(f+Sign(g,f))-h))/x
; следующее QR-преобразование
            c=1.0
            s=1.0
            For i1=l To k1
               i=i1+1
               g=rv1(i)
               y=w(i)
               h=s * g
               g=c * g
               z=Sqr(f * f+h * h)
               rv1(i1)=z
               c=f/z
               s=h/z
               f=x * c+g * s
               g=-x * s+g * c
               h=y * s
               y=y * c
               If matv
                  For j=1 To n
                     x=v(j,i1)
                     z=v(j,i)
                     v(j,i1)=x * c+z * s
                     v(j,i)=-x * s+z * c
                  Next
               EndIf
               z=Sqr(f * f+h * h)
               w(i1)=z
               If  z <> 0.0
                  c=f/z
                  s=h/z
               EndIf
               f=c * g+s * y
               x=-s * g+c * y
               If Not matu :    Continue : EndIf
               For j=1 To m
                  y=u(j,i1)
                  z=u(j,i)
                  u(j,i1)=y * c+z * s
                  u(j,i)=-y * s+z * c
               Next
            Next
            rv1(l)=0.0
            rv1(k)=f
            w(k)=x
         ForEver
; сходимость
         If z >= 0.0 :    Continue : EndIf
         w(k)=-z
         If Not matv :   Continue : EndIf
         For j=1 To n
            v(j,k)=-v(j,k)
         Next
       Next
       ProcedureReturn 1
EndProcedure

@ZAK53 · 13.11.2014, 02:53 **[ТС]**

Пропихнул по половинкам. чтобы прогнать тест, надо их поженить.

PureBasic
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
;------------------------------------------------------------------------------------------------------------------------
 
;              Регулярное решение линейной системы A * X=Y
;              на основе сингулярного разложения A=UWV'
 
 
;     SUBROUTINE REGSOLV(mn,m,n,U,W,V,Y,X,relerr)
;      integer*4 mn       ; строчный размер массивов U и V, объявленный в вызывающей программе !! На РВ этот параметр не нужен
;      integer*4 m        ; число строк матрицы U
;      integer*4 n        ; число столбцов матрицы U и порядок V
;      real*8    U(mn,n)  ; массив, содержащий матрицу U
;      real*8    W(n) ; массив, содержащий диагональ матрицы W
;      real*8    V(mn,n)  ; массив, содержащий матрицу V
;      real*8    Y(m)
;      real*8    X(n)     ; решение
;      real*8    relerr   ; граница относительной ошибки данных (Y).
 
Procedure RegSolv (m,n, Array U.d(2), Array W.d(1), Array V.d(2), Array Y.d(1), Array X.d(1), relerr.d)
 
      Protected.l i,j
      Protected.d rerr,tau,Wmax,rer,s
 
      rerr=1
      Repeat
        rerr=rerr/2
      Until 1+rerr=1
 
      rer=relerr
      If rer < rerr * n : rer=rerr * n  : EndIf
 
      Wmax=0.0
      For i=1 To n
         If W(i) > Wmax : Wmax=W(i)  : EndIf
         X(i)=0.0
      Next
 
      tau=rer * Wmax
 
      For j=1 To n
         If W(j) <= tau : Continue  : EndIf
         s=0.0
         For i=1 To m
            s=s+U(i,j) * Y(i)
         Next
         s=s/W(j)
         For i=1 To n
            X(i)=X(i)+s * V(i,j)
         Next
      Next
      ProcedureReturn
    EndProcedure
;---------------------------------------
Procedure.i max (a.i,b.i)
  Protected.i c
  If a => b : c=a : Else : c=b : EndIf
  ProcedureReturn c
EndProcedure
;-------------------------------------------------
 
;       TEST
 
Define.i i,j,nm, m=5, n=3, matu=1, matv=1, res, count
Define.d relerr=0
 
nm=max(m,n)
 
Dim U.d (nm,n)
Dim V.d (nm,n)
Dim W.d (n)
 
count=0
For i=1 To n
  For j=1 To m
    count=count+1
    U(j,i)=count
  Next
Next
Debug "Исходная матрица:"
For i=1 To m
  Debug StrD(U(i,1))+"  "+StrD(U(i,2))+"  "+StrD(U(i,3))
Next
 
res = Svd1(m, n, W(), matu, U(), matv, V())
 
Debug " "
Debug "SVD1 return " + Str(res)
Debug " "
Debug "Первые n столбцов матрицы U:"
For i=1 To m
  Debug StrD(U(i,1))+"  "+StrD(U(i,2))+"  "+StrD(U(i,3))
Next
Debug " "
Debug "Mатрица V:"
For i=1 To n
  Debug StrD(V(i,1))+"  "+StrD(V(i,2))+"  "+StrD(V(i,3))
Next
Debug " "
Debug "Диагональ матрицы W:"
Debug StrD(W(1))+"  "+StrD(W(2))+"  "+StrD(W(3))
 
Dim Y.d (m)
Dim X.d (n)
 
Debug " "
Debug "Положим свободный член системы уравнений Y=(5,5,5,5,5)'"
 
For i=1 To m
  Y(i)=5
Next
 
RegSolv (m,n, U(), W(), V(), Y(), X(), relerr)
 
Debug "Решение с минимальной нормой:"
Debug StrD(X(1))+"  "+StrD(X(2))+"  "+StrD(X(3))
MessageRequester ("X=", StrD(X(1))+"  "+StrD(X(2))+"  "+StrD(X(3)),0)
End

@gazlan · 13.11.2014, 03:15

До кучи. Оригинал.

Новые блоги и статьи Все статьи Все блоги /
SDL3 для Web (WebAssembly): Работа со звуком через SDL3_mixer 8Observer8 08.02.2026 Содержание блога Пошагово создадим проект для загрузки звукового файла и воспроизведения звука с помощью библиотеки SDL3_mixer. Звук будет воспроизводиться по клику мышки по холсту на Desktop и по. . .	SDL3 для Web (WebAssembly): Основы отладки веб-приложений на SDL3 по USB и Wi-Fi, запущенных в браузере мобильных устройств 8Observer8 07.02.2026 Содержание блога Браузер Chrome имеет средства для отладки мобильных веб-приложений по USB. В этой пошаговой инструкции ограничимся работой с консолью. Вывод в консоль - это часть процесса. . .	SDL3 для Web (WebAssembly): Обработчик клика мыши в браузере ПК и касания экрана в браузере на мобильном устройстве 8Observer8 02.02.2026 Содержание блога Для начала пошагово создадим рабочий пример для подготовки к экспериментам в браузере ПК и в браузере мобильного устройства. Потом напишем обработчик клика мыши и обработчик. . .	Философия технологии iceja 01.02.2026 На мой взгляд у человека в технических проектах остается роль генерального директора. Все остальное нейронки делают уже лучше человека. Они не могут нести предпринимательские риски, не могут. . .
SDL3 для Web (WebAssembly): Вывод текста со шрифтом TTF с помощью SDL3_ttf 8Observer8 01.02.2026 Содержание блога В этой пошаговой инструкции создадим с нуля веб-приложение, которое выводит текст в окне браузера. Запустим на Android на локальном сервере. Загрузим Release на бесплатный. . .	SDL3 для Web (WebAssembly): Сборка C/C++ проекта из консоли 8Observer8 30.01.2026 Содержание блога Если вы откроете примеры для начинающих на официальном репозитории SDL3 в папке: examples, то вы увидите, что все примеры используют следующие четыре обязательные функции, а. . .	SDL3 для Web (WebAssembly): Установка Emscripten SDK (emsdk) и CMake для сборки C и C++ приложений в Wasm 8Observer8 30.01.2026 Содержание блога Для того чтобы скачать Emscripten SDK (emsdk) необходимо сначало скачать и уставить Git: Install for Windows. Следуйте стандартной процедуре установки Git через установщик. . . .	SDL3 для Android: Подключение Box2D v3, физика и отрисовка коллайдеров 8Observer8 29.01.2026 Содержание блога Box2D - это библиотека для 2D физики для анимаций и игр. С её помощью можно определять были ли коллизии между конкретными объектами. Версия v3 была полностью переписана на Си, в. . .

@ZAK53 0 / 0 / 0 Регистрация: 11.11.2013 Сообщений: 17

	Сингулярное разложение матриц 13.11.2014, 00:26. Показов 3886. Ответов 3 Метки нет (Все метки) Не уверен, что это кому-нить нужно, но что делаю, о том пою. Сингулярное разложение (SVD) - очень полезный, а иногда и необходимый аппарат для считальщиков. Применить можно много для чего, в частности - для решения линейных систем с плохо обусловленными или вообще неполного ранга матрицами (тестовый пример как раз с такой матрицей). Еще - для решения переопределенных (уравнений больше, чем неизвестных) систем методом Наименьших Квадратов. В наш век, когда космические корабли бороздят и вычислительные мощности настольных ЭВМ огромны, вообще нет резона применять для решения систем более простые и экономные алгоритмы типа Гаусса. С ними можно получить "изумительные" результаты при плохой обусловленности. Ниже привожу перевод на РВ ФОРТРАН 4 -подпрограммы из книги [Форсайт, Малькольм, Моулер. Машинные методы математических вычислений. М., 1982(?)], которая, в свою очередь, являлась переводом с небольшой модификацией АЛГОЛ-процедуры [Уилкинсон, Райнш. Справочник алгоритмов на АЛГОЛ]. Тестовый пример тоже взят из [Форсайт, Малькольм, Моулер]. Описание оставил фортрановское ... Нулевые индексы в массивах не используются, индексация начинается с 1. Потамушто. Единственное, что порадовало в РВ, так это то ли отсутствие оптимизатора, то ли его ненавязчивость. Благодаря этому проходят тесты на пренебрежимую малость типа if A+anorm = anorm и оптимизатор не превращает их в if A=0. Черррт, исходник не проходит по размеру, попробую прицепить прицеп. Прицеп тоже не лезет, если кому интересно - возьмите с мого яндекс-диска: https://yadi.sk/d/wLnV4TYncg5iJ 0