| 150 | |
| 151 | @pytest.mark.skipif(nogo, reason="Missing modules (autograd or pymanopt)") |
| 152 | def test_ewca(): |
| 153 | d = 5 |
| 154 | n_samples = 50 |
| 155 | k = 3 |
| 156 | rng = np.random.RandomState(0) |
| 157 | |
| 158 | # generate gaussian dataset |
| 159 | A = rng.normal(size=(d, d)) |
| 160 | Q, _ = np.linalg.qr(A) |
| 161 | D = rng.normal(size=d) |
| 162 | D = (D / np.linalg.norm(D)) ** 4 |
| 163 | cov = Q @ np.diag(D) @ Q.T |
| 164 | X = rng.multivariate_normal(np.zeros(d), cov, size=n_samples) |
| 165 | X = X - X.mean(0, keepdims=True) |
| 166 | assert X.shape == (n_samples, d) |
| 167 | |
| 168 | # compute first 3 components with BCD |
| 169 | pi, U = ot.dr.ewca( |
| 170 | X, reg=0.01, method="BCD", k=k, verbose=1, sinkhorn_method="sinkhorn_log" |
| 171 | ) |
| 172 | assert pi.shape == (n_samples, n_samples) |
| 173 | assert (pi >= 0).all() |
| 174 | assert np.allclose(pi.sum(0), 1 / n_samples, atol=1e-3) |
| 175 | assert np.allclose(pi.sum(1), 1 / n_samples, atol=1e-3) |
| 176 | assert U.shape == (d, k) |
| 177 | assert np.allclose(U.T @ U, np.eye(k), atol=1e-3) |
| 178 | |
| 179 | # test that U contains the principal components |
| 180 | U_first_eigvec = np.linalg.svd(X.T, full_matrices=False)[0][:, :k] |
| 181 | _, cos, _ = np.linalg.svd(U.T @ U_first_eigvec, full_matrices=False) |
| 182 | assert np.allclose(cos, np.ones(k), atol=1e-3) |
| 183 | |
| 184 | # compute first 3 components with MM |
| 185 | pi, U = ot.dr.ewca( |
| 186 | X, reg=0.01, method="MM", k=k, verbose=1, sinkhorn_method="sinkhorn_log" |
| 187 | ) |
| 188 | assert pi.shape == (n_samples, n_samples) |
| 189 | assert (pi >= 0).all() |
| 190 | assert np.allclose(pi.sum(0), 1 / n_samples, atol=1e-3) |
| 191 | assert np.allclose(pi.sum(1), 1 / n_samples, atol=1e-3) |
| 192 | assert U.shape == (d, k) |
| 193 | assert np.allclose(U.T @ U, np.eye(k), atol=1e-3) |
| 194 | |
| 195 | # test that U contains the principal components |
| 196 | U_first_eigvec = np.linalg.svd(X.T, full_matrices=False)[0][:, :k] |
| 197 | _, cos, _ = np.linalg.svd(U.T @ U_first_eigvec, full_matrices=False) |
| 198 | assert np.allclose(cos, np.ones(k), atol=1e-3) |
| 199 | |
| 200 | # compute last 3 components |
| 201 | pi, U = ot.dr.ewca( |
| 202 | X, reg=100000, method="MM", k=k, verbose=1, sinkhorn_method="sinkhorn_log" |
| 203 | ) |
| 204 | |
| 205 | # test that U contains the last principal components |
| 206 | U_last_eigvec = np.linalg.svd(X.T, full_matrices=False)[0][:, -k:] |
| 207 | _, cos, _ = np.linalg.svd(U.T @ U_last_eigvec, full_matrices=False) |
| 208 | assert np.allclose(cos, np.ones(k), atol=1e-3) |