Related to this question: How to have negative zero always formatted as positive zero in a python string?
I have the following function that implements Matlab's orth.m using Numpy. I have a docstring test that relies on np.array2string using suppress_small=True, which will make small values round to zero. However, sometimes they round to positive zero and sometimes they round to negative zero, depending on whether the answer comes out as 1e-16 or -1e-17 or similar. Which case happens is based on the SVD decomposition, and can vary from platform to platform or across Python versions depending on which underlying linear algebra solver is used (BLAS, Lapack, etc.)
What's the best way to design the docstring test to account for this?
In the final doctest, sometimes the Q[0, 1] term is -0. and sometimes it's 0.
import doctest
import numpy as np
def orth(A):
r"""
Orthogonalization basis for the range of A.
That is, Q.T @ Q = I, the columns of Q span the same space as the columns of A, and the number
of columns of Q is the rank of A.
Parameters
----------
A : 2D ndarray
Input matrix
Returns
-------
Q : 2D ndarray
Orthogonalization matrix of A
Notes
-----
#. Based on the Matlab orth.m function.
Examples
--------
>>> import numpy as np
Full rank matrix
>>> A = np.array([[1, 0, 1], [-1, -2, 0], [0, 1, -1]])
>>> r = np.linalg.matrix_rank(A)
>>> print(r)
3
>>> Q = orth(A)
>>> with np.printoptions(precision=8):
... print(Q)
[[-0.12000026 -0.80971228 0.57442663]
[ 0.90175265 0.15312282 0.40422217]
[-0.41526149 0.5664975 0.71178541]]
Rank deficient matrix
>>> A = np.array([[1, 0, 1], [0, 1, 0], [1, 0, 1]])
>>> r = np.linalg.matrix_rank(A)
>>> print(r)
2
>>> Q = orth(A)
>>> print(np.array2string(Q, precision=8, suppress_small=True)) # Sometimes this fails
[[-0.70710678 -0. ]
[ 0. 1. ]
[-0.70710678 0. ]]
"""
# compute the SVD
(Q, S, _) = np.linalg.svd(A, full_matrices=False)
# calculate a tolerance based on the first eigenvalue (instead of just using a small number)
tol = np.max(A.shape) * S[0] * np.finfo(float).eps
# sum the number of eigenvalues that are greater than the calculated tolerance
r = np.sum(S > tol, axis=0)
# return the columns corresponding to the non-zero eigenvalues
Q = Q[:, np.arange(r)]
return Q
if __name__ == '__main__':
doctest.testmod(verbose=False)
You can print the rounded array plus 0.0 to eliminate the -0:
A = np.array([[1, 0, 1], [0, 1, 0], [1, 0, 1]])
Q = orth(A)
Q[0,1] = -1e-16 # simulate a small floating point deviation
print(np.array2string(Q.round(8)+0.0, precision=8, suppress_small=True))
#[[-0.70710678 0. ]
# [ 0. 1. ]
# [-0.70710678 0. ]]
So your doc string should be:
>>> Q = orth(A)
>>> print(np.array2string(Q.round(8)+0.0, precision=8, suppress_small=True)) # guarantee non-negative zeros
[[-0.70710678 0. ]
[ 0. 1. ]
[-0.70710678 0. ]]