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

@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

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

Новые блоги и статьи Все статьи Все блоги /
Отображение реквизитов в документе по условию и контроль их заполнения Maks 04.04.2026 Алгоритм из решения ниже реализован на примере нетипового документа "ПланированиеСпецтехники", разработанного в конфигурации КА2. Данный документ берёт данные из другого нетипового документа. . .	Фото всей Земли с борта корабля Orion миссии Artemis II kumehtar 04.04.2026 Это первое подобное фото сделанное человеком за 50 лет. Снимок называют новым вариантом легендарной фотографии «The Blue Marble» 1972 года, сделанной с борта корабля «Аполлон-17». Новое фото. . .	Вывод диалогового окна перед закрытием, если документ не проведён Maks 04.04.2026 Алгоритм из решения ниже реализован на примере нетипового документа "СписаниеМатериалов", разработанного в конфигурации КА2. Задача: реализовать программный контроль на предмет проведения документа. . .	Программный контроль заполнения реквизита табличной части документа Maks 02.04.2026 Алгоритм из решения ниже реализован на примере нетипового документа "СписаниеМатериалов", разработанного в конфигурации КА2. Задача: реализовать контроль заполнения реквизита "ПричинаСписания". . .
wmic не является внутренней или внешней командой Maks 02.04.2026 Решение: DISM / Online / Add-Capability / CapabilityName:WMIC~~~~ Отсюда: https:/ / winitpro. ru/ index. php/ 2025/ 02/ 14/ komanda-wmic-ne-naydena/	Программная установка даты и запрет ее изменения Maks 02.04.2026 Алгоритм из решения ниже реализован на примере нетипового документа "СписаниеМатериалов", разработанного в конфигурации КА2. Задача: при создании документов установить период списания автоматически. . .	Вывод данных в справочнике через динамический список Maks 01.04.2026 Реализация из решения ниже выполнена на примере нетипового справочника "Спецтехника" разработанного в конфигурации КА2. Задача: вывести данные из ТЧ нетипового документа. . .	Программное заполнения текстового поля в реквизите формы документа Maks 01.04.2026 Алгоритм из решения ниже реализован на нетиповом документе "ВыдачаОборудованияНаСпецтехнику" разработанного в конфигурации КА2, в дополнении к предыдущему решению. На форме документа создается. . .

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

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