1 use na::{DMatrix, Matrix6};
2
3 #[cfg(feature = "proptest-support")]
4 mod proptest_tests {
5 macro_rules! gen_tests(
6 ($module: ident, $scalar: expr, $scalar_type: ty) => {
7 mod $module {
8 use na::{
9 DMatrix, DVector, Matrix2, Matrix3, Matrix4,
10 ComplexField
11 };
12 use std::cmp;
13 #[allow(unused_imports)]
14 use crate::core::helper::{RandScalar, RandComplex};
15 use crate::proptest::*;
16 use proptest::{prop_assert, proptest};
17
18 proptest! {
19 #[test]
20 fn svd(m in dmatrix_($scalar)) {
21 let svd = m.clone().svd(true, true);
22 let recomp_m = svd.clone().recompose().unwrap();
23 let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
24 let ds = DMatrix::from_diagonal(&s.map(|e| ComplexField::from_real(e)));
25
26 prop_assert!(s.iter().all(|e| *e >= 0.0));
27 prop_assert!(relative_eq!(&u * ds * &v_t, recomp_m, epsilon = 1.0e-5));
28 prop_assert!(relative_eq!(m, recomp_m, epsilon = 1.0e-5));
29 }
30
31 #[test]
32 fn svd_static_5_3(m in matrix5x3_($scalar)) {
33 let svd = m.svd(true, true);
34 let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
35 let ds = Matrix3::from_diagonal(&s.map(|e| ComplexField::from_real(e)));
36
37 prop_assert!(s.iter().all(|e| *e >= 0.0));
38 prop_assert!(relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5));
39 prop_assert!(u.is_orthogonal(1.0e-5));
40 prop_assert!(v_t.is_orthogonal(1.0e-5));
41 }
42
43 #[test]
44 fn svd_static_5_2(m in matrix5x2_($scalar)) {
45 let svd = m.svd(true, true);
46 let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
47 let ds = Matrix2::from_diagonal(&s.map(|e| ComplexField::from_real(e)));
48
49 prop_assert!(s.iter().all(|e| *e >= 0.0));
50 prop_assert!(relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5));
51 prop_assert!(u.is_orthogonal(1.0e-5));
52 prop_assert!(v_t.is_orthogonal(1.0e-5));
53 }
54
55 #[test]
56 fn svd_static_3_5(m in matrix3x5_($scalar)) {
57 let svd = m.svd(true, true);
58 let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
59
60 let ds = Matrix3::from_diagonal(&s.map(|e| ComplexField::from_real(e)));
61
62 prop_assert!(s.iter().all(|e| *e >= 0.0));
63 prop_assert!(relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5));
64 }
65
66 #[test]
67 fn svd_static_2_5(m in matrix2x5_($scalar)) {
68 let svd = m.svd(true, true);
69 let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
70 let ds = Matrix2::from_diagonal(&s.map(|e| ComplexField::from_real(e)));
71
72 prop_assert!(s.iter().all(|e| *e >= 0.0));
73 prop_assert!(relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5));
74 }
75
76 #[test]
77 fn svd_static_square(m in matrix4_($scalar)) {
78 let svd = m.svd(true, true);
79 let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
80 let ds = Matrix4::from_diagonal(&s.map(|e| ComplexField::from_real(e)));
81
82 prop_assert!(s.iter().all(|e| *e >= 0.0));
83 prop_assert!(relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5));
84 prop_assert!(u.is_orthogonal(1.0e-5));
85 prop_assert!(v_t.is_orthogonal(1.0e-5));
86 }
87
88 #[test]
89 fn svd_static_square_2x2(m in matrix2_($scalar)) {
90 let svd = m.svd(true, true);
91 let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
92 let ds = Matrix2::from_diagonal(&s.map(|e| ComplexField::from_real(e)));
93
94 prop_assert!(s.iter().all(|e| *e >= 0.0));
95 prop_assert!(relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5));
96 prop_assert!(u.is_orthogonal(1.0e-5));
97 prop_assert!(v_t.is_orthogonal(1.0e-5));
98 }
99
100 #[test]
101 fn svd_pseudo_inverse(m in dmatrix_($scalar)) {
102 let svd = m.clone().svd(true, true);
103 let pinv = svd.pseudo_inverse(1.0e-10).unwrap();
104
105 if m.nrows() > m.ncols() {
106 prop_assert!((pinv * m).is_identity(1.0e-5))
107 } else {
108 prop_assert!((m * pinv).is_identity(1.0e-5))
109 }
110 }
111
112 #[test]
113 fn svd_solve(n in PROPTEST_MATRIX_DIM, nb in PROPTEST_MATRIX_DIM) {
114 let n = cmp::max(1, cmp::min(n, 10));
115 let nb = cmp::min(nb, 10);
116 let m = DMatrix::<$scalar_type>::new_random(n, n).map(|e| e.0);
117
118 let svd = m.clone().svd(true, true);
119
120 if svd.rank(1.0e-7) == n {
121 let b1 = DVector::<$scalar_type>::new_random(n).map(|e| e.0);
122 let b2 = DMatrix::<$scalar_type>::new_random(n, nb).map(|e| e.0);
123
124 let sol1 = svd.solve(&b1, 1.0e-7).unwrap();
125 let sol2 = svd.solve(&b2, 1.0e-7).unwrap();
126
127 let recomp = svd.recompose().unwrap();
128
129 prop_assert!(relative_eq!(m, recomp, epsilon = 1.0e-6));
130 prop_assert!(relative_eq!(&m * &sol1, b1, epsilon = 1.0e-6));
131 prop_assert!(relative_eq!(&m * &sol2, b2, epsilon = 1.0e-6));
132 }
133 }
134 }
135 }
136 }
137 );
138
139 gen_tests!(complex, complex_f64(), RandComplex<f64>);
140 gen_tests!(f64, PROPTEST_F64, RandScalar<f64>);
141 }
142
143 // Test proposed on the issue #176 of rulinalg.
144 #[test]
145 #[rustfmt::skip]
svd_singular()146 fn svd_singular() {
147 let m = DMatrix::from_row_slice(24, 24, &[
148 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
149 -1.0, -1.0, -1.0, -1.0, -1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
150 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
151 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0,
152 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 0.0, 1.0, 1.0, 1.0,
153 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
154 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
155 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
156 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
157 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
158 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
159 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
160 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
161 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
162 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
163 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
164 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
165 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
166 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 4.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
167 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
168 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
169 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
170 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0,
171 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0]);
172
173 let svd = m.clone().svd(true, true);
174 let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
175 let ds = DMatrix::from_diagonal(&s);
176
177 assert!(s.iter().all(|e| *e >= 0.0));
178 assert!(u.is_orthogonal(1.0e-5));
179 assert!(v_t.is_orthogonal(1.0e-5));
180 assert_relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5);
181 }
182
183 // Same as the previous test but with one additional row.
184 #[test]
185 #[rustfmt::skip]
svd_singular_vertical()186 fn svd_singular_vertical() {
187 let m = DMatrix::from_row_slice(25, 24, &[
188 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
189 -1.0, -1.0, -1.0, -1.0, -1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
190 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
191 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0,
192 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 0.0, 1.0, 1.0, 1.0,
193 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
194 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
195 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
196 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
197 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
198 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
199 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
200 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
201 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
202 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
203 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
204 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
205 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
206 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 4.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
207 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
208 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
209 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
210 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0,
211 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0,
212 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0]);
213
214
215 let svd = m.clone().svd(true, true);
216 let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
217 let ds = DMatrix::from_diagonal(&s);
218
219 assert!(s.iter().all(|e| *e >= 0.0));
220 assert_relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5);
221 }
222
223 // Same as the previous test but with one additional column.
224 #[test]
225 #[rustfmt::skip]
svd_singular_horizontal()226 fn svd_singular_horizontal() {
227 let m = DMatrix::from_row_slice(24, 25, &[
228 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
229 -1.0, -1.0, -1.0, -1.0, -1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
230 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
231 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
232 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 0.0, 1.0, 1.0, 1.0, 1.0,
233 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
234 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
235 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
236 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
237 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
238 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
239 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
240 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
241 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
242 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
243 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
244 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
245 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
246 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 4.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
247 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
248 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
249 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
250 0.0, 0.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
251 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, -4.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]);
252
253 let svd = m.clone().svd(true, true);
254 let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
255 let ds = DMatrix::from_diagonal(&s);
256
257 assert!(s.iter().all(|e| *e >= 0.0));
258 assert_relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5);
259 }
260
261 #[test]
svd_zeros()262 fn svd_zeros() {
263 let m = DMatrix::from_element(10, 10, 0.0);
264 let svd = m.clone().svd(true, true);
265 assert_eq!(Ok(m), svd.recompose());
266 }
267
268 #[test]
svd_identity()269 fn svd_identity() {
270 let m = DMatrix::<f64>::identity(10, 10);
271 let svd = m.clone().svd(true, true);
272 assert_eq!(Ok(m), svd.recompose());
273
274 let m = DMatrix::<f64>::identity(10, 15);
275 let svd = m.clone().svd(true, true);
276 assert_eq!(Ok(m), svd.recompose());
277
278 let m = DMatrix::<f64>::identity(15, 10);
279 let svd = m.clone().svd(true, true);
280 assert_eq!(Ok(m), svd.recompose());
281 }
282
283 #[test]
284 #[rustfmt::skip]
svd_with_delimited_subproblem()285 fn svd_with_delimited_subproblem() {
286 let mut m = DMatrix::<f64>::from_element(10, 10, 0.0);
287 m[(0, 0)] = 1.0; m[(0, 1)] = 2.0;
288 m[(1, 1)] = 0.0; m[(1, 2)] = 3.0;
289 m[(2, 2)] = 4.0; m[(2, 3)] = 5.0;
290 m[(3, 3)] = 6.0; m[(3, 4)] = 0.0;
291 m[(4, 4)] = 8.0; m[(3, 5)] = 9.0;
292 m[(5, 5)] = 10.0; m[(3, 6)] = 11.0;
293 m[(6, 6)] = 12.0; m[(3, 7)] = 12.0;
294 m[(7, 7)] = 14.0; m[(3, 8)] = 13.0;
295 m[(8, 8)] = 16.0; m[(3, 9)] = 17.0;
296 m[(9, 9)] = 18.0;
297 let svd = m.clone().svd(true, true);
298 assert_relative_eq!(m, svd.recompose().unwrap(), epsilon = 1.0e-7);
299
300 // Rectangular versions.
301 let mut m = DMatrix::<f64>::from_element(15, 10, 0.0);
302 m[(0, 0)] = 1.0; m[(0, 1)] = 2.0;
303 m[(1, 1)] = 0.0; m[(1, 2)] = 3.0;
304 m[(2, 2)] = 4.0; m[(2, 3)] = 5.0;
305 m[(3, 3)] = 6.0; m[(3, 4)] = 0.0;
306 m[(4, 4)] = 8.0; m[(3, 5)] = 9.0;
307 m[(5, 5)] = 10.0; m[(3, 6)] = 11.0;
308 m[(6, 6)] = 12.0; m[(3, 7)] = 12.0;
309 m[(7, 7)] = 14.0; m[(3, 8)] = 13.0;
310 m[(8, 8)] = 16.0; m[(3, 9)] = 17.0;
311 m[(9, 9)] = 18.0;
312 let svd = m.clone().svd(true, true);
313 assert_relative_eq!(m, svd.recompose().unwrap(), epsilon = 1.0e-7);
314
315 let svd = m.transpose().svd(true, true);
316 assert_relative_eq!(m.transpose(), svd.recompose().unwrap(), epsilon = 1.0e-7);
317 }
318
319 #[test]
320 #[rustfmt::skip]
svd_fail()321 fn svd_fail() {
322 let m = Matrix6::new(
323 0.9299319121545955, 0.9955870335651049, 0.8824725266413644, 0.28966880207132295, 0.06102723649846409, 0.9311880746048009,
324 0.5938395242304351, 0.8398522876024204, 0.06672831951963198, 0.9941213119963099, 0.9431846038057834, 0.8159885168706427,
325 0.9121962883152357, 0.6471119669367571, 0.4823309702814407, 0.6420516076705516, 0.7731203925207113, 0.7424069470756647,
326 0.07311092531259344, 0.5579247949052946, 0.14518764691585773, 0.03502980663114896, 0.7991329455957719, 0.4929930019965745,
327 0.12293810556077789, 0.6617084679545999, 0.9002240700227326, 0.027153062135304884, 0.3630189466989524, 0.18207502727558866,
328 0.843196731466686, 0.08951878746549924, 0.7533450877576973, 0.009558876499740077, 0.9429679490873482, 0.9355764454129878);
329 let svd = m.clone().svd(true, true);
330 let recomp = svd.recompose().unwrap();
331 assert_relative_eq!(m, recomp, epsilon = 1.0e-5);
332 }
333
334 #[test]
svd_err()335 fn svd_err() {
336 let m = DMatrix::from_element(10, 10, 0.0);
337 let svd = m.clone().svd(false, false);
338 assert_eq!(
339 Err("SVD recomposition: U and V^t have not been computed."),
340 svd.clone().recompose()
341 );
342 assert_eq!(
343 Err("SVD pseudo inverse: the epsilon must be non-negative."),
344 svd.clone().pseudo_inverse(-1.0)
345 );
346 }
347