"""Utility module of reusable algebra routines."""
import math
import geomstats.backend as gs
EPSILON = 1e-6
COS_TAYLOR_COEFFS = [
1.0,
-1.0 / math.factorial(2),
+1.0 / math.factorial(4),
-1.0 / math.factorial(6),
+1.0 / math.factorial(8),
]
SINC_TAYLOR_COEFFS = [
1.0,
-1.0 / math.factorial(3),
+1.0 / math.factorial(5),
-1.0 / math.factorial(7),
+1.0 / math.factorial(9),
]
INV_SINC_TAYLOR_COEFFS = [1, 1.0 / 6.0, 7.0 / 360.0, 31.0 / 15120.0, 127.0 / 604800.0]
INV_TANC_TAYLOR_COEFFS = [1.0, -1.0 / 3.0, -1.0 / 45.0, -2.0 / 945.0, -1.0 / 4725.0]
COSC_TAYLOR_COEFFS = [
1.0 / 2.0,
-1.0 / math.factorial(4),
+1.0 / math.factorial(6),
-1.0 / math.factorial(8),
+1.0 / math.factorial(10),
]
VAR_INV_TAN_TAYLOR_COEFFS = [1.0 / 12.0, 1.0 / 720.0, 1.0 / 30240.0, 1.0 / 1209600.0]
SINHC_TAYLOR_COEFFS = [
1.0,
1 / math.factorial(3),
1 / math.factorial(5),
1 / math.factorial(7),
1 / math.factorial(9),
]
COSH_TAYLOR_COEFFS = [
1.0,
1 / math.factorial(2),
1 / math.factorial(4),
1 / math.factorial(6),
1 / math.factorial(8),
]
INV_SINHC_TAYLOR_COEFFS = [
1.0,
-1.0 / 6.0,
7.0 / 360.0,
-31.0 / 15120.0,
127.0 / 604800.0,
]
INV_TANH_TAYLOR_COEFFS = [1.0, 1.0 / 3.0, -1.0 / 45.0, 2.0 / 945.0, -1.0 / 4725.0]
ARCTANH_CARD_TAYLOR_COEFFS = [1.0, 1.0 / 3.0, 1.0 / 5.0, 1 / 7.0, 1.0 / 9]
cos_close_0 = {"function": gs.cos, "coefficients": COS_TAYLOR_COEFFS}
sinc_close_0 = {"function": lambda x: gs.sin(x) / x, "coefficients": SINC_TAYLOR_COEFFS}
inv_sinc_close_0 = {
"function": lambda x: x / gs.sin(x),
"coefficients": INV_SINC_TAYLOR_COEFFS,
}
inv_tanc_close_0 = {
"function": lambda x: x / gs.tan(x),
"coefficients": INV_TANC_TAYLOR_COEFFS,
}
cosc_close_0 = {
"function": lambda x: (1 - gs.cos(x)) / x**2,
"coefficients": COSC_TAYLOR_COEFFS,
}
var_sinc_close_0 = {
"function": lambda x: (x - gs.sin(x)) / x**3,
"coefficients": [-k for k in SINC_TAYLOR_COEFFS[1:]],
}
var_inv_tanc_close_0 = {
"function": lambda x: (1 - (x / gs.tan(x))) / x**2,
"coefficients": VAR_INV_TAN_TAYLOR_COEFFS,
}
sinch_close_0 = {
"function": lambda x: gs.sinh(x) / x,
"coefficients": SINHC_TAYLOR_COEFFS,
}
cosh_close_0 = {"function": gs.cosh, "coefficients": COSH_TAYLOR_COEFFS}
inv_sinch_close_0 = {
"function": lambda x: x / gs.sinh(x),
"coefficients": INV_SINHC_TAYLOR_COEFFS,
}
inv_tanh_close_0 = {
"function": lambda x: x / gs.tanh(x),
"coefficients": INV_TANH_TAYLOR_COEFFS,
}
arctanh_card_close_0 = {
"function": lambda x: gs.arctanh(x) / x,
"coefficients": ARCTANH_CARD_TAYLOR_COEFFS,
}
[docs]
def from_vector_to_diagonal_matrix(vector, num_diag=0):
"""Create diagonal matrices from rows of a matrix.
Parameters
----------
vector : array-like, shape=[m, n]
num_diag : int
number of diagonal in result matrix. If 0, the result matrix is a
diagonal matrix; if positive, the result matrix has an upper-right
non-zero diagonal; if negative, the result matrix has a lower-left
non-zero diagonal.
Optional, Default: 0.
Returns
-------
diagonals : array-like, shape=[m, n, n]
3-dimensional array where the `i`-th n-by-n array `diagonals[i, :, :]`
is a diagonal matrix containing the `i`-th row of `vector`.
"""
num_columns = gs.shape(vector)[-1]
identity = gs.eye(num_columns, dtype=vector.dtype)
diagonals = gs.einsum("...i,ij->...ij", vector, identity)
diagonals = gs.to_ndarray(diagonals, to_ndim=3)
num_lines = diagonals.shape[0]
if num_diag > 0:
left_zeros = gs.zeros((num_lines, num_columns, num_diag))
lower_zeros = gs.zeros((num_lines, num_diag, num_columns + num_diag))
diagonals = gs.concatenate((left_zeros, diagonals), axis=2)
diagonals = gs.concatenate((diagonals, lower_zeros), axis=1)
elif num_diag < 0:
num_diag = gs.abs(num_diag)
right_zeros = gs.zeros((num_lines, num_columns, num_diag))
upper_zeros = gs.zeros((num_lines, num_diag, num_columns + num_diag))
diagonals = gs.concatenate((diagonals, right_zeros), axis=2)
diagonals = gs.concatenate((upper_zeros, diagonals), axis=1)
return gs.squeeze(diagonals) if gs.ndim(vector) == 1 else diagonals
[docs]
def taylor_exp_even_func(point, taylor_function, order=5, tol=EPSILON):
"""Taylor Approximation of an even function around zero.
Parameters
----------
point : array-like
Argument of the function to approximate.
taylor_function : dict with following keys
function : callable
Even function to approximate around zero.
coefficients : list
Taylor coefficients of even order at zero.
order : int
Order of the Taylor approximation.
Optional, Default: 5.
tol : float
Threshold to use the approximation instead of the function's value.
Where `abs(point) <= tol`, the approximation is returned.
Returns
-------
function_value: array-like
Value of the function at point.
"""
approx = gs.einsum(
"k,k...->...",
gs.array(taylor_function["coefficients"][:order]),
gs.array([point**k for k in range(order)]),
)
point_ = gs.where(gs.abs(point) <= tol, tol, point)
exact = taylor_function["function"](gs.sqrt(point_))
result = gs.where(gs.abs(point) < tol, approx, exact)
return result
[docs]
def flip_determinant(matrix, det):
"""Change sign of the determinant if it is negative.
For a batch of matrices, multiply the matrices which have negative
determinant by a diagonal matrix :math:`diag(1,...,1,-1) from the right.
This changes the sign of the last column of the matrix.
Parameters
----------
matrix : array-like, shape=[...,n ,m]
Matrix to transform.
det : array-like, shape=[...]
Determinant of matrix, or any other scalar to use as threshold to
determine whether to change the sign of the last column of matrix.
Returns
-------
matrix_flipped : array-like, shape=[..., n, m]
Matrix with the sign of last column changed if det < 0.
"""
if gs.any(det < 0):
ones = gs.ones(matrix.shape[-1])
reflection_vec = gs.concatenate([ones[:-1], gs.array([-1.0])], axis=0)
mask = gs.cast(det < 0, matrix.dtype)
sign = mask[..., None] * reflection_vec + (1.0 - mask)[..., None] * ones
return gs.einsum("...ij,...j->...ij", matrix, sign)
return matrix
[docs]
def rotate_points(points, end_point):
"""Apply to points the rotation from north_pole to end_point.
A QR decomposition is used to find the rotation that maps the north pole
(1, 0,...,0) to the end_point, then this rotation is applied to the
input points.
Parameters
----------
points : array-like, shape=[..., n]
Points to rotate.
end_point : array-like, shape=[n, ]
Point to parametrise the rotation.
Returns
-------
rotated_points : array-like, shape=[..., n]
Points after the rotation.
"""
n = end_point.shape[0]
base_point = gs.array([1.0] + [0] * (n - 1))
embedded = gs.concatenate([end_point[None, :], gs.zeros((n - 1, n))])
norm = gs.linalg.norm(end_point)
q, _ = gs.linalg.qr(gs.transpose(embedded) / norm)
new_points = gs.matmul(points[None, :], gs.transpose(q)) * norm
if not gs.allclose(gs.matmul(q, base_point[:, None])[:, 0], end_point):
new_points = -new_points
return new_points[0]
[docs]
def columnwise_scaling(vec, mat):
r"""Column-wise scaling.
Equivalent to :math:`AD`, where :math:`D` is a
diagonal matrix.
Parameters
----------
vec : array-like, shape=[..., k]
Vector of scalings.
mat :array-like, shape=[..., n, k]
Matrix.
Returns
-------
column_scaled_mat : array-like, shape=[..., n, k]
"""
return vec[..., None, :] * mat
[docs]
def rowwise_scaling(vec, mat):
r"""Row-wise scaling.
Equivalent to :math:`DA`, where :math:`D` is a
diagonal matrix.
Parameters
----------
vec : array-like, shape=[..., n]
Vector of scalings.
mat :array-like, shape=[..., n, k]
Matrix.
Returns
-------
row_scaled_mat : array-like, shape=[..., n, k]
"""
return vec[..., :, None] * mat