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
|
#! /usr/bin/env python
from __future__ import print_function
from openturns import *
TESTPREAMBLE()
def quadM(m, n):
res = Matrix(m, n)
for i in range(m):
for j in range(n):
res[i, j] = (i + 1.0) ** (j + 1.0)
return res
def testQR(m, n, full, keep):
matrix1 = quadM(m, n)
print("M=", matrix1)
Q, R = matrix1.computeQR(full, keep)
print("full=", full, "keep=", keep)
print('Q= ', Q)
print('R=', R)
print('Q*R=', Q * R)
if keep:
print('M2=', matrix1)
try:
# Square case
matrix1 = quadM(3, 3)
matrix1.setName("matrix1")
print("matrix1 = ", repr(matrix1))
result1 = matrix1.computeSingularValues()
print("svd (svd only)= ", repr(result1))
result1, u, v = matrix1.computeSVD(True)
print("svd (svd + U, V full)= ", repr(result1))
result1, u, v = matrix1.computeSVD(False)
print("svd (svd + U, V small)= ", repr(result1),
", U=", repr(u), ", v=", repr(v))
# Rectangular case, m < n
matrix1 = quadM(3, 5)
matrix1.setName("matrix1")
print("matrix1 = ", repr(matrix1))
result1 = matrix1.computeSingularValues()
print("svd (svd only)= ", repr(result1))
result1, u, v = matrix1.computeSVD(True)
print("svd (svd + U, V full)= ", repr(result1))
result1, u, v = matrix1.computeSVD(False)
print("svd (svd + U, V small)= ", repr(result1),
", U=", repr(u), ", v=", repr(v))
# Rectangular case, m > n
matrix1 = quadM(5, 3)
matrix1.setName("matrix1")
print("matrix1 = ", repr(matrix1))
result1 = matrix1.computeSingularValues()
print("svd (svd only)= ", repr(result1))
result1, u, v = matrix1.computeSVD(True)
print("svd (svd + U, V full)= ", repr(result1))
#result1, u, v = matrix1.computeSVD(False)
# print "svd (svd + U, V small)= ", repr(result1), ", U=", repr(u), ",
# v=", repr(v)
for iFull in range(2):
for iKeep in range(2):
testQR(3, 3, iFull == 1, iKeep == 1)
testQR(3, 5, iFull == 1, iKeep == 1)
testQR(5, 3, iFull == 1, iKeep == 1)
except:
import sys
print("t_Matrix_decomposition.py", sys.exc_info()[0], sys.exc_info()[1])
|
