import numpy as np
from ddg.math.linalg import rank
def _rnd_matrix_full_rank(m, n, seed=None, atol=1e-13, rtol=None):
"""Helper function for random_matrix"""
gen = np.random.default_rng(seed=seed)
A = gen.standard_normal(size=(m, n))
r = rank(A, atol=atol, rtol=rtol)
while r < min(m, n):
A = gen.standard_normal(size=(m, n))
r = rank(A, atol=atol, rtol=rtol)
return A
[docs]def random_matrix(m, n, rank=None, seed=None, atol=1e-13, rtol=None):
"""Returns a random (m x n)-matrix of a specific rank.
Parameters
----------
m : int
Number of rows of matrix.
n : int
Number of columns of matrix.
rank : int, optional.
Rank of matrix. If None, a full rank matrix
seed : int, numpy.random.Generator or None (default=None)
Passes a fixed seed to default Generator class. If passed a Generator,
just uses that generator for random generation.
atol : float (default=1e-13)
rtol : float (default=None)
The absolute and relative tolerance used in the computation of rank of the
matrix.
This function uses the global tolerance defaults if `atol` or `rtol`
are set to None. See :py:mod:`ddg.nonexact` for details.
Raises
------
ValueError
If rank of matrix is strictly greater than min(m, n)
Returns
-------
numpy.ndarray
A random (m x n)-matrix
Notes
-----
"""
if rank and rank != min(m, n):
if rank > min(m, n):
raise ValueError(
"Can not return matrix with rank greater than "
"the number of its rows or columns."
)
else:
A = _rnd_matrix_full_rank(m, rank, seed=seed, atol=atol, rtol=rtol)
gen = np.random.default_rng(seed=seed)
B = np.array(
[
gen.standard_normal() * A[:, gen.integers(0, A.shape[1])]
for _ in range(0, n - rank)
]
).T
return np.hstack((A, B))
else:
return _rnd_matrix_full_rank(m, n, seed=seed, atol=atol, rtol=rtol)
[docs]def random_matrix_close_to_identity(n, seed=None):
"""Returns a random (n x n)-matrix close to identity matrix.
Parameters
----------
n : int
Dimension of the square matrix.
seed : int, numpy.random.Generator or None (default=None)
Passes a fixed seed to default Generator class. If passed a Generator,
just uses that generator for random generation.
Returns
-------
numpy.ndarray
A (n x n)-matrix
"""
gen = np.random.default_rng(seed=seed)
A = gen.standard_normal(size=(n, n))
Id = np.identity(n)
B = Id + A / (n**2 * np.linalg.norm(A))
return B
[docs]def random_point_on_sphere(n, seed=None):
"""Draws a point on the n-dimensional unit sphere with
uniform distribution.
Parameters
----------
n : int
Dimension of the unit sphere.
seed : float, optional.
Passes a fixed seed to default Generator class.
Returns
-------
numpy.array of shape (n+1,)
Point on the n-dimensional unit sphere.
"""
generator = np.random.default_rng(seed=seed)
coordinates = generator.normal(0.0, 1.0, size=(n + 1,))
return coordinates / np.linalg.norm(coordinates)