hypothesis_helpers.invertible_matrices generates numerically singular matrices

invertible_matrices is:

array-api-tests/array_api_tests/hypothesis_helpers.py

Lines 426 to 443 in ef5b39f

    
           def invertible_matrices(draw, dtypes=floating_dtypes, stack_shapes=shapes()): 
        
               # For now, just generate stacks of diagonal matrices. 
        
               stack_shape = draw(stack_shapes) 
        
               n = draw(integers(0, SQRT_MAX_ARRAY_SIZE // max(math.prod(stack_shape), 1)),) 
        
               dtype = draw(dtypes) 
        
               elements = one_of( 
        
                   from_dtype(dtype, min_value=0.5, allow_nan=False, allow_infinity=False), 
        
                   from_dtype(dtype, max_value=-0.5, allow_nan=False, allow_infinity=False), 
        
               ) 
        
               d = draw(arrays(dtype, shape=(*stack_shape, 1, n), elements=elements)) 
        
               # Functions that require invertible matrices may do anything when it is 
        
               # singular, including raising an exception, so we make sure the diagonals 
        
               # are sufficiently nonzero to avoid any numerical issues. 
        
               assert xp.all(xp.abs(d) >= 0.5) 
        
               diag_mask = xp.arange(n) == xp.reshape(xp.arange(n), (n, 1)) 
        
               return xp.where(diag_mask, d, xp.zeros_like(d))

There is no upper bound on the magnitude of the entries, so it is free to generate (and reliably generates, in my experience) very poorly conditioned matrices that mpmath considers to be numerically singular.

The comment says that it is trying to avoid numerical issues; to do that, it would need to bound the relative magnitudes of the entries.

The assertion xp.all(xp.abs(d) >= 0.5) can also be removed since that is guaranteed by the parameters provided to one_of.

Similarly, the allowed powers of n when testing matrix_power can result in numerical difficulties, depending on the algorithm.

array-api-tests/array_api_tests/test_linalg.py

Lines 481 to 504 in ef5b39f

    
           matrix_power_n = shared(integers(-100, 100), key='matrix_power n') 
        
           @pytest.mark.unvectorized 
        
           @pytest.mark.xp_extension('linalg') 
        
           @given( 
        
               # Generate any square matrix if n >= 0 but only invertible matrices if n < 0 
        
               x=matrix_power_n.flatmap(lambda n: invertible_matrices() if n < 0 else 
        
                                        arrays(dtype=all_floating_dtypes(), 
        
                                                   shape=square_matrix_shapes)), 
        
               n=matrix_power_n, 
        
           ) 
        
           def test_matrix_power(x, n): 
        
               res = linalg.matrix_power(x, n) 
        
               ph.assert_dtype("matrix_power", in_dtype=x.dtype, out_dtype=res.dtype) 
        
               ph.assert_result_shape("matrix_power", in_shapes=[x.shape], 
        
                                      out_shape=res.shape, expected=x.shape) 
        
               if n == 0: 
        
                   true_val = lambda x: xp.eye(x.shape[0], dtype=x.dtype) 
        
               else: 
        
                   true_val = None 
        
               # _test_stacks only works with array arguments 
        
               func = lambda x: linalg.matrix_power(x, n) 
        
               _test_stacks(func, x, res=res, true_val=true_val)

Roughly, I think the choice of n should ensure that the condition number of the expected result is much less than 1/eps.

Yeah, linalg testing barely does any value testing, and even trying to do more of it is either very messy or outright impossible at the array-api-tests level. Some prior discussions: https://github.com/data-apis/array-api-tests/pull/314/changes#r1859369692

Yeah, I've noticed. But this is just about the input. The test does run the function, and the numerically singular matrix causes an error, which causes test failure.

Ouch. What array library, how did you catch the error? And would you be interested in sending a PR with a fix?

Surprisingly MPArray. Not that important, I know, but maybe it will save some other library frustration down the line. Yeah, I can submit a PR.

Navigation Menu

Repository navigation

`hypothesis_helpers.invertible_matrices` generates numerically singular matrices #467

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

	def invertible_matrices(draw, dtypes=floating_dtypes, stack_shapes=shapes()):
	# For now, just generate stacks of diagonal matrices.
	stack_shape = draw(stack_shapes)
	n = draw(integers(0, SQRT_MAX_ARRAY_SIZE // max(math.prod(stack_shape), 1)),)
	dtype = draw(dtypes)
	elements = one_of(
	from_dtype(dtype, min_value=0.5, allow_nan=False, allow_infinity=False),
	from_dtype(dtype, max_value=-0.5, allow_nan=False, allow_infinity=False),
	)
	d = draw(arrays(dtype, shape=(*stack_shape, 1, n), elements=elements))

	# Functions that require invertible matrices may do anything when it is
	# singular, including raising an exception, so we make sure the diagonals
	# are sufficiently nonzero to avoid any numerical issues.
	assert xp.all(xp.abs(d) >= 0.5)

	diag_mask = xp.arange(n) == xp.reshape(xp.arange(n), (n, 1))
	return xp.where(diag_mask, d, xp.zeros_like(d))

	matrix_power_n = shared(integers(-100, 100), key='matrix_power n')
	@pytest.mark.unvectorized
	@pytest.mark.xp_extension('linalg')
	@given(
	# Generate any square matrix if n >= 0 but only invertible matrices if n < 0
	x=matrix_power_n.flatmap(lambda n: invertible_matrices() if n < 0 else
	arrays(dtype=all_floating_dtypes(),
	shape=square_matrix_shapes)),
	n=matrix_power_n,
	)
	def test_matrix_power(x, n):
	res = linalg.matrix_power(x, n)

	ph.assert_dtype("matrix_power", in_dtype=x.dtype, out_dtype=res.dtype)
	ph.assert_result_shape("matrix_power", in_shapes=[x.shape],
	out_shape=res.shape, expected=x.shape)

	if n == 0:
	true_val = lambda x: xp.eye(x.shape[0], dtype=x.dtype)
	else:
	true_val = None
	# _test_stacks only works with array arguments
	func = lambda x: linalg.matrix_power(x, n)
	_test_stacks(func, x, res=res, true_val=true_val)

hypothesis_helpers.invertible_matrices generates numerically singular matrices #467

Description

Activity

mdhaber commented on Sep 14, 2026

ev-br commented on Sep 14, 2026

mdhaber commented on Sep 14, 2026

ev-br commented on Sep 14, 2026

mdhaber commented on Sep 14, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions

`hypothesis_helpers.invertible_matrices` generates numerically singular matrices #467