BaiduFMX

1
unit SimpleSVD;
2

3
interface
4

5
uses
6
  System.Classes,
7
  System.SysUtils,
8
  System.Math;
9

10
type
11
  TSVDFloat = Double;
12
  TSVDFloatArray = TArray<TSVDFloat>;
13

14
  TMatrix = class
15
  private
16
    FRows: Integer;
17
    FCols: Integer;
18
    FName: string;
19
    FType: Integer;
20
    FMat: TSVDFloatArray;
21
    function GetMat(Index: Integer): TSVDFloat;
22
    procedure SetMat(Index: Integer; const Value: TSVDFloat);
23
    function GetMatrix(row, col: Integer): TSVDFloat;
24
    procedure SetMatrix(row, col: Integer; const Value: TSVDFloat);
25
  public const
26
    AS_MATRIX = 1;
27
    AS_VECTOR = 2;
28
  public
29
    constructor Create(M, N: Integer; aName: string; aType: Integer);
30
    procedure Output(Lines: TStrings);
31
    procedure SetAll(const aMat: TSVDFloatArray);
32
    property Mat[Index: Integer]: TSVDFloat read GetMat write SetMat;
33
    property Rows: Integer read FRows;
34
    property Cols: Integer read FCols;
35
    property Name: string read FName;
36
    property MatrixType: Integer read FType;
37
    property Matrix[row, col: Integer]: TSVDFloat read GetMatrix write SetMatrix; default;
38
  end;
39

40
  TSVD = class
41
  public
42
    constructor Create;
43
    function Multiply (m1, m2, y: TMatrix): Integer;
44
    function PinvCompute(A: TMatrix; rows, cols: Integer): TMatrix;
45
    function PseudoInverseSVD(const V,S,U: TSVDFloatArray; rows,cols: Integer): TSVDFloatArray;
46
    procedure GivensL(var S_: TSVDFloatArray; const dim: TArray<Integer>; m: Integer; a, b: TSVDFloat);
47
    procedure GivensR(var S_: TSVDFloatArray; const dim: TArray<Integer>; m: Integer; a, b: TSVDFloat);
48
    procedure SVDDecomp(const dim: TArray<Integer>; var U_, S_, V_: TSVDFloatArray; eps: Double);
49
    procedure SVDCompute(M, N: Integer; const A: TSVDFloatArray; LDA: Integer;
50
      var S, U: TSVDFloatArray; LDU: Integer; var VT: TSVDFloatArray;
51
      LDVT: Integer);
52
  end;
53
implementation
54

55
{ TMatrix }
56

57
constructor TMatrix.Create(M, N: Integer; aName: string; aType: Integer);
58
begin
59
  FRows := M;
60
  FCols := N;
61
  FName := aName;
62
  FType := aType;
63
  SetLength(FMat, M * N);
64
end;
65

66
function TMatrix.GetMat(Index: Integer): TSVDFloat;
67
begin
68
  Assert(Index < FRows * FCols);
69
  Result := FMat[Index];
70
end;
71

72
function TMatrix.GetMatrix(row, col: Integer): TSVDFloat;
73
begin
74
  Assert((row < FRows) and (col < FCols));
75
  Result := FMat[row * FCols + col];
76
end;
77

78
procedure TMatrix.Output(Lines: TStrings);
79
var
80
  s: string;
81
  i, j: Integer;
82
begin
83
  if(FType=AS_MATRIX) then
84
  begin
85
    Lines.Add('Matrix ' + FName + ':');
86
    for i:=0 to FRows-1 do
87
    begin
88
      s := '';
89
      for j:=0 to FCols-1 do
90
      begin
91
        if(FMat[i*cols+j]<10E-13) and (FMat[i*cols+j]>-10E-13) then
92
            mat[i*cols+j] := 0;
93
        s := s + Format('%12.4g', [FMat[i*FCols+j]]);
94
      end;
95
      Lines.Add(s);
96
    end;
97
    Lines.Add('');
98
  end
99
  else if(FType=AS_VECTOR) then
100
  begin
101
    Lines.Add('Vector ' + FName + ':');
102
    s := '';
103
    for i:=0 to FRows-1 do
104
    begin
105
      if(FMat[i]<10E-13) and (FMat[i]>-10E-13) then
106
        mat[i*cols] := 0;
107
      s := s + Format('%13.4f', [FMat[i]]);
108
    end;
109
    Lines.Add(s);
110
    Lines.Add('');
111
  end;
112
end;
113

114
procedure TMatrix.SetAll(const aMat: TSVDFloatArray);
115
begin
116
  FMat := aMat;
117
end;
118

119
procedure TMatrix.SetMat(Index: Integer; const Value: TSVDFloat);
120
begin
121
  Assert(Index < FRows * FCols);
122
  FMat[Index] := Value;
123
end;
124

125
procedure TMatrix.SetMatrix(row, col: Integer; const Value: TSVDFloat);
126
begin
127
  Assert((row < FRows) and (col < FCols));
128
  FMat[row * FCols + col] := Value;
129
end;
130

131
{ TSVD }
132

133
constructor TSVD.Create;
134
begin
135

136
end;
137

138
//#define U(i,j) U_[(i)*dim[0]+(j)]
139
//#define S(i,j) S_[(i)*dim[1]+(j)]
140
//#define V(i,j) V_[(i)*dim[1]+(j)]
141

142
procedure TSVD.GivensL(var S_: TSVDFloatArray; const dim: TArray<Integer>;
143
  m: Integer; a, b: TSVDFloat);
144
var
145
  r, c, s: TSVDFloat;
146
  i, idx0, idx1: Integer;
147
  S0, S1: TSVDFloat;
148
begin
149
  r:=sqrt(a*a+b*b);
150
  c:=a/r;
151
  s:=-b/r;
152

153
  for i:=0 to dim[1]-1 do
154
  begin
155
    idx0 := (m+0)*dim[0]+(i);  //S(m+0,i)
156
    idx1 := (m+1)*dim[0]+(i);  //S(m+1,i)
157
    S0:=S_[idx0];
158
    S1:=S_[idx1];
159
    S_[idx0] := S_[idx0] + S0*(c-1) + S1*(-s);
160
    S_[idx1] := S_[idx1] + S0*(s) + S1*(c-1);
161
  end;
162
end;
163

164
procedure TSVD.GivensR(var S_: TSVDFloatArray; const dim: TArray<Integer>;
165
  m: Integer; a, b: TSVDFloat);
166
var
167
  r, c, s: TSVDFloat;
168
  i, idx0, idx1: Integer;
169
  S0, S1: TSVDFloat;
170
begin
171
  r:=sqrt(a*a+b*b);
172
  c:=a/r;
173
  s:=-b/r;
174
  for i:=0 to dim[0]-1 do
175
  begin
176
    idx0 := (i)*dim[0]+(m+0);  //S(i,m+0)
177
    idx1 := (i)*dim[0]+(m+1);  //S(i,m+1);
178

179
    S0:=S_[idx0];
180
    S1:=S_[idx1];
181
    S_[idx0] := S_[idx0] + S0*(c-1) + S1*(-s);
182
    S_[idx1] := S_[idx1] + S0*(s) + S1*(c-1);
183
  end;
184
end;
185

186
function TSVD.Multiply(m1, m2, y: TMatrix): Integer;
187
var
188
  sum: Double;
189
  I, J, K: Integer;
190
begin
191
  if (m1.Cols <> m2.Rows) then
192
    Exit(0);
193

194
  for I := 0 to m1.Rows-1 do
195
  begin
196
    for J := 0 to m2.Cols-1 do
197
    begin
198
      sum := 0;
199
      for K := 0 to m1.Cols-1 do
200
      begin
201
        sum := sum + m1[i,k] * m2[k,j];
202
      end;
203
      y[i,j] := sum;
204
    end;
205
  end;
206

207
  Result := 1;
208
end;
209

210
function TSVD.PinvCompute(A: TMatrix; rows, cols: Integer): TMatrix;
211
var
212
  U, VT, S: TMatrix;
213
  Mat: TSVDFloatArray;
214
  I, J: Integer;
215
  m, n: Integer;
216
  tU, tS, tVT, pinv: TSVDFloatArray;
217
begin
218
  Result := TMatrix.Create(cols,rows,'A_pinv',TMatrix.AS_MATRIX);
219
  U:=TMatrix.Create(rows,cols,'U',TMatrix.AS_MATRIX);
220
  VT:=TMatrix.Create(cols,cols,'VT',TMatrix.AS_MATRIX);
221
  S:=TMatrix.Create(cols,cols,'S',TMatrix.AS_MATRIX);
222
  SetLength(Mat, rows*cols);
223
  try
224
    for I:=0 to rows*cols-1 do
225
      Mat[i]:=A.Mat[i];
226

227
    m := cols;
228
    n := rows;
229
    SetLength(tU, m*m);
230
    SetLength(tS, m);
231
    SetLength(tVT, m*n);
232

233
    SvdCompute(m, n, Mat, m, tS, tU, m, tVT, m);
234

235
    for I:=0 to rows-1 do
236
      for j:=0 to cols-1 do
237
        U.Mat[i*cols+j] := tVT[i*cols+j];
238

239
    for i:=0 to cols-1 do
240
      for j:=0 to cols-1 do
241
        VT.Mat[i*cols+j] := tU[i*cols+j];
242

243
    for i:=0 to cols-1 do
244
      for j:=0 to cols-1 do
245
      begin
246
        if(i=j) then
247
          S.Mat[i*cols+j]:=tS[i]
248
        else
249
          S.Mat[i*cols+j]:=0;
250
      end;
251
    pinv := PseudoInverseSVD(tVT,tS,tU,m,n);
252
    Result.SetAll(pinv);
253
  finally
254
    S.Free;
255
    VT.Free;
256
    U.Free;
257
  end;
258
end;
259

260
function TSVD.PseudoInverseSVD(const V, S, U: TSVDFloatArray; rows,
261
  cols: Integer): TSVDFloatArray;
262
var
263
  temp: TSVDFloatArray;
264
  I, J, K: Integer;
265
  sum: TSVDFloat;
266
begin
267
  SetLength(temp, rows*rows);
268

269
  for I:=0 to rows-1 do
270
    for J := 0 to rows-1 do
271
      temp[i*rows+j] := U[j*rows+i]/S[j];
272
  SetLength(Result, rows*cols);
273

274
  for i:=0 to rows-1 do
275
    for j:=0 to cols-1 do
276
    begin
277
      sum:=0;
278
      for k:=0 to rows-1 do
279
        sum:=sum+temp[i*rows+k]*V[j*rows+k];
280
      Result[i*cols+j]:=sum;
281
    end;
282
end;
283

284
procedure TSVD.SVDCompute(M, N: Integer; const A: TSVDFloatArray; LDA: Integer;
285
  var S, U: TSVDFloatArray; LDU: Integer; var VT: TSVDFloatArray;
286
  LDVT: Integer);
287
var
288
  dim: TArray<Integer>;
289
  U_, V_, S_: TSVDFloatArray;
290
  i, j: Integer;
291
  tmp: TSVDFloat;
292
begin
293
  dim := TArray<Integer>.Create(Max(N, M), Min(N, M));
294
  SetLength(U_, dim[0]*dim[0]);
295
  SetLength(V_, dim[1]*dim[1]);
296
  SetLength(S_, dim[0]*dim[1]);
297
  FillChar(U_[0], dim[0]*dim[0]*SizeOf(TSVDFloat), 0);
298
  FillChar(V_[0], dim[1]*dim[1]*SizeOf(TSVDFloat), 0);
299

300
  if(dim[1]=M) then
301
  begin
302
    for I := 0 to dim[0]-1 do
303
      for J := 0 to dim[1]-1 do
304
      begin
305
        S_[i*dim[1]+j]:=A[i*lda+j];
306
      end;
307
  end
308
  else
309
  begin
310
    for I := 0 to dim[0]-1 do
311
      for J := 0 to dim[1]-1 do
312
      begin
313
        S_[i*dim[1]+j]:=A[j*lda+i];
314
      end;
315
  end;
316
  for i:=0 to dim[0]-1 do
317
  begin
318
    U_[i*dim[0]+i]:=1;
319
  end;
320
  for i:=0 to dim[1]-1 do
321
  begin
322
    V_[i*dim[1]+i]:=1;
323
  end;
324

325
  SvdDecomp(dim, U_, S_, V_, -1.0);
326

327
  for i:=0 to dim[1]-1 do
328
  begin
329
    S[i]:=S_[i*dim[1]+i];
330
  end;
331
  if(dim[1]=M) then
332
  begin
333
    for I := 0 to dim[1]-1 do
334
      for J := 0 to M-1 do
335
      begin
336
        if S[i] < 0.0 then
337
          tmp := -1.0
338
        else
339
          tmp := 1.0;
340
        U[j+ldu*i]:=V_[j+i*dim[1]]*tmp;
341
      end;
342
  end
343
  else
344
  begin
345
   for I := 0 to dim[1]-1 do
346
      for J := 0 to M-1 do
347
      begin
348
        if S[i] < 0.0 then
349
          tmp := -1.0
350
        else
351
          tmp := 1.0;
352
        U[j+ldu*i]:=U_[i+j*dim[0]]*tmp;
353
      end;
354
  end;
355
  if(dim[0]=N) then
356
  begin
357
    for i:=0 to N-1 do
358
      for j:=0 to dim[1]-1 do
359
      begin
360
        VT[j+ldvt*i]:=U_[j+i*dim[0]];
361
      end;
362
  end
363
  else
364
  begin
365
    for i:=0 to N-1 do
366
      for j:=0 to dim[1]-1 do
367
      begin
368
        VT[j+ldvt*i]:=V_[i+j*dim[1]];
369
      end;
370
  end;
371
  for i:=0 to dim[1]-1 do
372
  begin
373
    if S[i] < 0.0 then
374
      tmp := -1.0
375
    else
376
      tmp := 1.0;
377
    S[i]:=S[i]*tmp;
378
  end;
379
end;
380

381
procedure TSVD.SVDDecomp(const dim: TArray<Integer>; var U_, S_,
382
  V_: TSVDFloatArray; eps: Double);
383
var
384
  N: Integer;
385
  x1: TSVDFloat;
386
  x_inv_norm, dot_prod: Double;
387
  house_vec: TArray<Double>;
388
  i, j, K, idx, k0, idx2, i0, i1: Integer;
389
  tmp: TSVDFloat;
390
  alpha, beta, S_max: Double;
391
  c: array [0 .. 1, 0 .. 1] of Double;
392
  b, c2, d, lambda1, lambda2, d1, d2, mu: Double;
393
  dimU, dimV: TArray<Integer>;
394
begin
395
  Assert(dim[0] >= dim[1]);
396
  N := Min(dim[0], dim[1]);
397

398
  SetLength(house_vec, Max(dim[0], dim[1]));
399

400
  // S(i,j) S_[(i)*dim[1]+(j)]
401
  for i := 0 to N - 1 do
402
  begin
403
    x1 := S_[(i) * dim[1] + (i)]; // S(i,i);
404
    if (x1 < 0) then
405
      x1 := -x1;
406

407
    x_inv_norm := 0;
408

409
    for j := i to dim[0] - 1 do
410
    begin
411
      x_inv_norm := x_inv_norm + Sqr(S_[(j) * dim[1] + (i)]); // s(j,i)*s(j,i)
412
    end;
413
    if (x_inv_norm > 0) then
414
      x_inv_norm := 1 / sqrt(x_inv_norm);
415

416
    alpha := sqrt(1 + x1 * x_inv_norm);
417
    beta := x_inv_norm / alpha;
418

419
    house_vec[i] := -alpha;
420
    for j := i + 1 to dim[0] - 1 do
421
    begin
422
      house_vec[j] := -beta * S_[(j) * dim[1] + (i)];
423
    end;
424

425
    if (S_[(i) * dim[1] + (i)] < 0) then
426
      for j := i + 1 to dim[0] - 1 do
427
      begin
428
        house_vec[j] := -house_vec[j];
429
      end;
430

431
    for K := i to dim[1] - 1 do
432
    begin
433
      dot_prod := 0;
434
      for j := i to dim[0] - 1 do
435
      begin
436
        dot_prod := dot_prod + S_[(j) * dim[1] + (K)] * house_vec[j];
437
      end;
438
      for j := i to dim[0] - 1 do
439
      begin
440
        idx := (j) * dim[1] + (K);
441
        S_[idx] := S_[idx] - dot_prod * house_vec[j];
442
      end;
443
    end;
444

445
    // #define U(i,j) U_[(i)*dim[0]+(j)]
446
    for K := 0 to dim[0] - 1 do
447
    begin
448
      dot_prod := 0;
449
      for j := i to dim[0] - 1 do
450
      begin
451
        dot_prod := dot_prod + U_[(K) * dim[0] + (j)] * house_vec[j];
452
      end;
453
      for j := i to dim[0] - 1 do
454
      begin
455
        idx := (K) * dim[0] + (j);
456
        U_[idx] := U_[idx] - dot_prod * house_vec[j];
457
      end;
458
    end;
459

460
    if (i >= N - 1) then
461
      continue;
462
    begin
463
      x1 := S_[(i) * dim[1] + (i + 1)]; // S(i,i+1);
464
      if (x1 < 0) then
465
        x1 := -x1;
466

467
      x_inv_norm := 0;
468
      for j := i + 1 to dim[1] - 1 do
469
      begin
470
        x_inv_norm := x_inv_norm + Sqr(S_[(i) * dim[1] + (j)]); // S(i,j);
471
      end;
472
      if (x_inv_norm > 0) then
473
        x_inv_norm := 1 / sqrt(x_inv_norm);
474

475
      alpha := sqrt(1 + x1 * x_inv_norm);
476
      beta := x_inv_norm / alpha;
477

478
      house_vec[i + 1] := -alpha;
479
      for j := i + 2 to dim[1] - 1 do
480
      begin
481
        house_vec[j] := -beta * S_[(i) * dim[1] + (j)]; // S(i,j);
482
      end;
483
      if (S_[(i) * dim[1] + (i + 1)] < 0) then
484
        for j := i + 2 to dim[1] - 1 do
485
        begin
486
          house_vec[j] := -house_vec[j];
487
        end;
488
    end;
489

490
    for K := i to dim[0] - 1 do
491
    begin
492
      dot_prod := 0;
493
      for j := i + 1 to dim[1] - 1 do
494
      begin
495
        dot_prod := dot_prod + S_[(K) * dim[1] + (j)] * house_vec[j];
496
      end;
497
      for j := i + 1 to dim[1] - 1 do
498
      begin
499
        idx := (K) * dim[1] + (j);
500
        S_[idx] := S_[idx] - dot_prod * house_vec[j];
501
      end;
502
    end;
503

504
    // V(i,j) V_[(i)*dim[1]+(j)]
505
    for K := 0 to dim[1] - 1 do
506
    begin
507
      dot_prod := 0;
508
      for j := i + 1 to dim[1] - 1 do
509
      begin
510
        idx := (j) * dim[1] + (K);
511
        dot_prod := dot_prod + V_[idx] * house_vec[j];
512
      end;
513
      for j := i + 1 to dim[1] - 1 do
514
      begin
515
        idx := (j) * dim[1] + (K);
516
        V_[idx] := V_[idx] - dot_prod * house_vec[j];
517
      end;
518
    end;
519
  end;
520

521
  k0 := 0;
522
  if (eps < 0) then
523
  begin
524
    eps := 1.0;
525
    while (eps + 1.0 > 1.0) do
526
      eps := eps * 0.5;
527
    eps := eps * 64.0;
528
  end;
529

530
  while (k0 < dim[1] - 1) do
531
  begin
532
    S_max := 0.0;
533
    for i := 0 to dim[1] - 1 do
534
    begin
535
      idx := (i) * dim[1] + (i);
536
      S_max := Max(S_max, S_[idx]);
537
    end;
538

539
    while (k0 < dim[1] - 1) and (Abs(S_[(k0) * dim[1] + (k0 + 1)]) <= eps * S_max) do
540
      Inc(k0);
541
    if (k0 = (dim[1] - 1)) then
542
      continue;
543

544
    N := k0 + 2;
545
    while (N < dim[1]) and (Abs(S_[(N - 1) * dim[1] + (N)]) > eps * S_max) do
546
      Inc(N);
547

548
    begin
549
      idx := (N - 2) * dim[1] + (N - 2);
550
      c[0][0] := Sqr(S_[idx]);
551
      if (N - k0 > 2) then
552
      begin
553
        idx := (N - 3) * dim[1] + (N - 2);
554
        c[0][0] := c[0][0] + Sqr(S_[idx]);
555
      end;
556
      idx := (N - 2) * dim[1] + (N - 2);
557
      idx2 := (N - 2) * dim[1] + (N - 1);
558
      c[0][1] := S_[idx] * S_[idx2];
559
      c[1][0] := c[0][1];
560

561
      idx := (N - 1) * dim[1] + (N - 1);
562
      c[1][1] := Sqr(S_[idx]) + Sqr(S_[idx2]);
563

564
      b := -(c[0][0] + c[1][1]) / 2;
565
      c2 := c[0][0] * c[1][1] - c[0][1] * c[1][0];
566
      d := 0;
567
      if (b * b - c2 > 0) then
568
        d := sqrt(b * b - c2)
569
      else
570
      begin
571
        b := (c[0][0] - c[1][1]) / 2;
572
        c2 := -c[0][1] * c[1][0];
573
        if (b * b - c2 > 0) then
574
          d := sqrt(b * b - c2);
575
      end;
576

577
      lambda1 := -b + d;
578
      lambda2 := -b - d;
579

580
      d1 := lambda1 - c[1][1];
581
      if d1 < 0 then
582
        d1 := -d1;
583
      d2 := lambda2 - c[1][1];
584
      if d2 < 0 then
585
        d2 := -d2;
586
      if d1 < d2 then
587
        mu := lambda1
588
      else
589
        mu := lambda2;
590

591
      idx := (k0) * dim[1] + (k0);
592
      idx2 := (k0) * dim[1] + (k0 + 1);
593
      alpha := Sqr(S_[idx]) - mu;
594
      beta := S_[idx] * S_[idx2];
595
    end;
596

597
    for K := k0 to N - 2 do
598
    begin
599
      dimU := TArray<Integer>.Create(dim[0], dim[0]);
600
      dimV := TArray<Integer>.Create(dim[1], dim[1]);
601
      GivensR(S_, dim, K, alpha, beta);
602
      GivensL(V_, dimV, K, alpha, beta);
603

604
      idx := (K) * dim[1] + (K);
605
      idx2 := (K + 1) * dim[1] + (K);
606
      alpha := S_[idx];
607
      beta := S_[idx2];
608
      GivensL(S_, dim, K, alpha, beta);
609
      GivensR(U_, dimU, K, alpha, beta);
610

611
      idx := (K) * dim[1] + (K + 1);
612
      alpha := S_[idx];
613
      idx := (K) * dim[1] + (K + 2);
614
      beta := S_[idx];
615
    end;
616

617
    begin
618
      for i0 := k0 to N - 2 do
619
      begin
620
        for i1 := 0 to dim[1] - 1 do
621
        begin
622
          if (i0 > i1) or (i0 + 1 < i1) then
623
          begin
624
            idx := (i0) * dim[1] + (i1);
625
            S_[idx] := 0;
626
          end;
627
        end;
628
      end;
629
      for i0 := 0 to dim[0] - 1 do
630
      begin
631
        for i1 := k0 to N - 2 do
632
        begin
633
          if (i0 > i1) or (i0 + 1 < i1) then
634
          begin
635
            idx := (i0) * dim[1] + (i1);
636
            S_[idx] := 0;
637
          end;
638
        end;
639
      end;
640
      for i := 0 to dim[1] - 2 do
641
      begin
642
        idx := (i) * dim[1] + (i + 1);
643
        if (Abs(S_[idx]) <= eps * S_max) then
644
        begin
645
          S_[idx] := 0;
646
        end;
647
      end;
648
    end;
649
  end;
650
end;
651

652
end.
653