Skip to content

SVD decomposition float32 float64 stability #13309

@massich

Description

@massich

While working with #13303 we had issues with the randomized_svd stability.

So we start with the one from scipy and we found that this is not passing

    n_samples = 100
    n_features = 500
    rank = 5
    k = 10

    decimal = 5
    dtype = np.dtype(np.float64)

    # generate a matrix X of approximate effective rank `rank` and no noise
    # component (very structured signal):
    X = make_low_rank_matrix(n_samples=n_samples, n_features=n_features,
                             effective_rank=rank, tail_strength=0.0,
                             random_state=0).astype(dtype, copy=False)
    assert_equal(X.shape, (n_samples, n_features))

    # compute the singular values of X using the slow exact method
    U_64, s_64, V_64 = linalg.svd(X.astype(np.float64), full_matrices=False)
    U_32, s_32, V_32 = linalg.svd(X.astype(np.float32), full_matrices=False)

    # Assert types
    for current_var in (U_64, s_64, V_64):
        assert current_var.dtype == np.float64

    for current_var in (U_32, s_32, V_32):
        assert current_var.dtype == np.float32

    # Assert consistency
    assert_allclose(U_32, U_64, rtol=1e-5)
    assert_allclose(s_32, s_64, rtol=1e-5)
    assert_allclose(V_32, V_64, rtol=1e-5)

cc: @glemaitre, @thibsej

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions