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
|
# Author:: Ernest Ellingson
# Copyright:: Copyright (c) 2005
# These are extensions to the std-lib 'matrix' to allow an all ruby SVD
require 'matrix'
require 'mathn'
class Array
def sum(identity = 0, &block)
return identity unless size > 0
if block_given?
map(&block).sum
else
reduce(:+)
end
end
end
class Vector
def magnitude
sumsqs = 0.0
self.size.times do |i|
sumsqs += self[i] ** 2.0
end
Math.sqrt(sumsqs)
end
def normalize
nv = []
mag = self.magnitude
self.size.times do |i|
nv << (self[i] / mag)
end
Vector[*nv]
end
end
class Matrix
def Matrix.diag(s)
Matrix.diagonal(*s)
end
alias :trans :transpose
def SV_decomp(maxSweeps = 20)
if self.row_size >= self.column_size
q = self.trans * self
else
q = self * self.trans
end
qrot = q.dup
v = Matrix.identity(q.row_size)
azrot = nil
mzrot = nil
cnt = 0
s_old = nil
mu = nil
while true do
cnt += 1
for row in (0...qrot.row_size-1) do
for col in (1..qrot.row_size-1) do
next if row == col
h = Math.atan((2 * qrot[row,col])/(qrot[row,row]-qrot[col,col]))/2.0
hcos = Math.cos(h)
hsin = Math.sin(h)
mzrot = Matrix.identity(qrot.row_size)
mzrot[row,row] = hcos
mzrot[row,col] = -hsin
mzrot[col,row] = hsin
mzrot[col,col] = hcos
qrot = mzrot.trans * qrot * mzrot
v = v * mzrot
end
end
s_old = qrot.dup if cnt == 1
sum_qrot = 0.0
if cnt > 1
qrot.row_size.times do |r|
sum_qrot += (qrot[r,r]-s_old[r,r]).abs if (qrot[r,r]-s_old[r,r]).abs > 0.001
end
s_old = qrot.dup
end
break if (sum_qrot <= 0.001 and cnt > 1) or cnt >= maxSweeps
end # of do while true
s = []
qrot.row_size.times do |r|
s << Math.sqrt(qrot[r,r])
end
#puts "cnt = #{cnt}"
if self.row_size >= self.column_size
mu = self * v * Matrix.diagonal(*s).inverse
return [mu, v, s]
else
puts v.row_size
puts v.column_size
puts self.row_size
puts self.column_size
puts s.size
mu = (self.trans * v * Matrix.diagonal(*s).inverse)
return [mu, v, s]
end
end
def []=(i,j,val)
@rows[i][j] = val
end
end
|
