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
|
# File 'lib/classifier/extensions/vector.rb', line 52
def SV_decomp(max_sweeps = 20)
q_matrix = if row_size >= column_size
trans * self
else
self * trans
end
q_rotation_matrix = q_matrix.dup
v_matrix = Matrix.identity(q_matrix.row_size)
iteration_count = 0
previous_s_matrix = nil
loop do
iteration_count += 1
(0...(q_rotation_matrix.row_size - 1)).each do |row|
(1..(q_rotation_matrix.row_size - 1)).each do |col|
next if row == col
numerator = 2.0 * q_rotation_matrix[row, col]
denominator = q_rotation_matrix[row, row] - q_rotation_matrix[col, col]
angle = if denominator.abs < Vector::EPSILON
numerator >= 0 ? Math::PI / 4.0 : -Math::PI / 4.0
else
Math.atan(numerator / denominator) / 2.0
end
cosine = Math.cos(angle)
sine = Math.sin(angle)
rotation_matrix = Matrix.identity(q_rotation_matrix.row_size)
rotation_matrix[row, row] = cosine
rotation_matrix[row, col] = -sine
rotation_matrix[col, row] = sine
rotation_matrix[col, col] = cosine
q_rotation_matrix = rotation_matrix.trans * q_rotation_matrix * rotation_matrix
v_matrix *= rotation_matrix
end
end
previous_s_matrix = q_rotation_matrix.dup if iteration_count == 1
sum_of_differences = 0.to_r
if iteration_count > 1
q_rotation_matrix.row_size.times do |r|
difference = (q_rotation_matrix[r, r] - previous_s_matrix[r, r]).abs
sum_of_differences += difference.to_r if difference > 0.001
end
previous_s_matrix = q_rotation_matrix.dup
end
break if (sum_of_differences <= 0.001 && iteration_count > 1) || iteration_count >= max_sweeps
end
singular_values = q_rotation_matrix.row_size.times.map do |r|
Math.sqrt([q_rotation_matrix[r, r].to_f, 0.0].max)
end
safe_singular_values = singular_values.map { |v| [v, Vector::EPSILON].max }
u_matrix = (row_size >= column_size ? self : trans) * v_matrix * Matrix.diagonal(*safe_singular_values).inverse
[u_matrix, v_matrix, singular_values]
end
|