Generating covariance matrices (and rotations and affine transformations)

[JohannesBuchner]

Motivation

For scientific python (and perhaps computer vision as well), a very common application is linear algebra. The simplest objects of interest there are vectors, which are easy to generate following the hypothesis documentation.

The second most common object of interest are perhaps rotation matrices and covariance matrices. Covariance matrices are positive semi-definite matrices, so simply generating a matrix and then checking whether it is valid is not an efficient strategy. Searching "hypothesis covariance matrix" brings up literature on a quite different topic.

I had this problem recently, so I thought I would share a strategy to generate covariance matrices.

Strategy

The strategy is not surprising to those familiar with linear algebra:

1. generate eigenvalues and eigenvectors. Eigenvectors need to be orthogonal, which can be achieved with the Gram-Schmidt process (QR-factorization).
2. build the covariance matrix (eigvec @ eigval @ eigvec.T). A affine transformation matrix instead could be built with T = eigvec * eigval**-0.5. The q from QR factorizatoin are a rotation matrix.

@st.composite
def mean_and_cov(draw):
    dim = draw(st.integers(min_value=1, max_value=10))  # Arbitrary dimensionality
    mu = draw(arrays(np.float64, (dim,), elements=st.floats(-10, 10)))  # Mean vector
    eigval = draw(arrays(np.float64, (dim,), elements=st.floats(1e-6, 10)))  # Eigenvalues
    vectors = draw(arrays(np.float64, (dim,dim), elements=st.floats(-10, 10)).filter(valid_QR))  # Eigenvectros
    cov = make_covariance_matrix_via_QR(eigval, vectors)
    return dim, mu, cov

def make_covariance_matrix_via_QR(normalisations, vectors):
    q, r = np.linalg.qr(vectors)
    orthogonal_vectors = q @ np.diag(np.diag(r))
    cov = orthogonal_vectors @ np.diag(normalisations) @ orthogonal_vectors.T
    return cov

3. Nevertheless, for numerical reasons, this can still rarely produce matrices that cannot be inverted. So finally, in the test using it I verify that the matrix is valid.

def valid_covariance_matrix(A, min_std):
    if not np.isfinite(A).all():
        return False
    if (np.diag(A) <= min_std).any():
        return False

    try:
        np.linalg.inv(A)
    except np.linalg.LinAlgError:
        return False
    try:
        scipy.stats.multivariate_normal(mean=np.zeros(len(A)), cov=A)
    except ValueError:
        return False
    return True

@given(mean_and_cov())
def test_single(mean_cov):
    ndim, mu, cov = mean_cov
    if not valid_covariance_matrix(cov, min_std=1e-6):
        return
    assert mu.shape == (ndim,), (mu, mu.shape, ndim)
    assert cov.shape == (ndim,ndim), (cov, cov.shape, ndim)

Limitations

I am a beginner in hypothesis, so probably this can be written much better. For example, I don't understand how I can chain a strategy, generating ndim first, and then passing that into mean_and_cov?
In mean_and_cov, there are range constraints hard-coded that may need to be adjusted depending on the application. These could perhaps be parameters of the strategy.

Proposal

Perhaps this can be incorporated into hypothesis.extra.numpy.

Alternatives

An noteworthy alternative is to generate covariance matrices from a Wishart distribution:

    seed = draw(st.integers(min_value=1, max_value=100000))
    cov = scipy.stats.wishart.rvs(df=dim, scale=np.eye(dim), random_state=np.random.RandomState(seed)).reshape((dim, dim))

However, the drawback here is that hypothesis will have a hard time shrinking to similar simpler examples.

[Zac-HD]

Nice!

This is clearly very useful for people working in the relevant domains; I'm just not sure whether Hypothesis itself is the best place to put it, or whether a third-party extension might be better for maintainence (hypothesis-linalg? or as part of another scientific computing package, ala xarray.testing?). That's mostly because it seems unlikely that covariance matrices are the only such special arrays that we'd want to generate, but adding all of them to the general hypothesis.extra.numpy (or ...arrays) namespace would get quite crowded.

Random technical notes:

• Check out the implementation of the arrays() strategy to see how the shapes and elements arguments are implemented. Also note however that accepting "value or strategy" is against our API style guide, and only allowed for arrays() for backwards-compatibility reasons.
• I'd recommend implementing this against the array-api strategies rather than numpy strategies, for flexibility
• You can apply valid_covariance_matrix() as a filter, or better yet as an assume() call inside the mean_and_cov() strategy. The @st.composite decorator, or .flatmap() method, make it easy to 'chain' strategies together.
• Shrinking is a pretty important feature for most users, and so it's usually worth going to a fair bit more trouble to make this work well. It's not only useful when shrinking a failing test either; the same principles which make the correspondence between underlying choices and high-level data and behavior also make trying out variations more effective in e.g. coverage-guided fuzzing.

[JohannesBuchner]

Thank you for the comment @Zac-HD.

I am afraid I am overloaded already by maintaining some dozen projects on pypi, so I don't think I can start a third-party extension project at this point. I am not sure I understood how to implement the random technical notes yet, it will take me some time. But I wanted to say thank you for taking the time and care to respond.

[JohannesBuchner]

FYI, this lead to scikit-learn/scikit-learn#29989 and https://github.com/JohannesBuchner/gmm-tests