import numpy as np
from ddg.abc import NonExact
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]@NonExact.nonexact_function
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
is returned
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
The absolute tolerance used in the computation
of rank of the matrix.
rtol : float
The relative tolerance used in the computation
of rank of the matrix.
Raises
------
ValueError
If rank of matrix is strictly greater than min(m, n)
Returns
-------
numpy.ndarray
A random (m x n)-matrix
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol`
are set to None. See ddg.abc.NonExact for details.
"""
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.
"""
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.,1.,size=(n+1,))
return coordinates/np.linalg.norm(coordinates)