import numpy as np
from typing import Optional, Union, List
from itertools import combinations
from .base import Vector, Matrix
from .tolerances import tolerances
[docs]
def range_space(A: Matrix):
r"""
Compute the range-space of a matrix.
"""
# Perform SVD
U, S, Vt = np.linalg.svd(A, full_matrices=False)
# Compute the number of non-zero singular values
rank = np.sum(S > np.finfo(S.dtype).eps)
# The range space is spanned by the first 'rank' columns of U
range_space = U[:, :rank]
return range_space
# def detect_dependent_rows(matrix: Matrix, zero_tol: Optional[float] = tolerances['zero']) -> IntVector:
# """
# Detect the row indices of linear dependent rows in a matrix.
# """
# _, S, _ = np.linalg.svd(matrix, full_matrices=False)
# dependent_indices = np.where(S <= zero_tol)[0]
# return dependent_indices
# def remove_dependent_rows(matrix: Matrix, zero_tol: Optional[float] = tolerances['zero']) -> Matrix:
# """
# Remove linearly dependent rows in a matrix.
# """
# _, S, _ = np.linalg.svd(matrix, full_matrices=False)
# independent_indices = np.where(S > zero_tol)[0]
# independent_rows = matrix[independent_indices, :]
# return independent_rows
[docs]
def is_lexico_positive(matrix: Union[Vector, Matrix], zero_tol: float = tolerances['zero']) -> bool:
r""" Tell if a matrix or vector is lexico positive.
- A vector is lexico positive if its first non-zero elements is nonnegative (following the style in MPT3).
- A matrix is lexico positive if all its rows are lexico positive.
Args:
matrix: a ndarray of size [m,n]. The array can be a 1d array representing a vector.
zero_tol: tolerance for a number to be recognized as non-zero.
Returns:
a bool variable implying the matrix is lexico positive or not.
"""
# process the dimension of inputs
if matrix.size == 0:
raise ValueError("The input matrix is empty!")
if len(matrix.shape) > 2:
raise ValueError("Wrong array size!")
if len(matrix.shape) == 1:
matrix = np.reshape(matrix, (1, matrix.size))
if zero_tol <= 0:
raise ValueError("zero_tol must be a strictly positive number")
pos = np.cumsum(matrix > zero_tol, axis=1).flatten()
neg = np.cumsum(matrix < -zero_tol, axis=1).flatten()
return not np.any((pos - neg) < 0)
[docs]
def lexico_argmax(matrix: Matrix, zero_tol: float = tolerances['zero']) -> List[int]:
r"""Find the (row) index of the lexico maximum of a matrix.
The lexico maximum of a set of vectors :math:`\{a_1, a_2, ..., a_m\}` is the vector :math:`a_j` satisfying
that :math:`a_j - a_i` is lexico positive for all :math:`i \in \{1,2,...,m\}/\{j\}`.
Args:
matrix: ndarray.
zero_tol: the tolerance to determine if two real number are identical or not.
Returns:
An integer vector consisting of index of lexico maximum row (the index could be non-unique).
"""
return lexico_argmin(-matrix, zero_tol)
[docs]
def lexico_argmin(matrix: Matrix, zero_tol: float = tolerances['zero']) -> List[int]:
r"""Find the (row) index of the lexico minimum of a matrix.
The lexico minimum of a set of vectors :math:`\{a_1, a_2, ..., a_m\}` is the vector :math:`a_j` satisfying
that :math:`a_i - a_j` is lexico positive for all :math:`i \in \{1,2,...,m\}/\{j\}`.
Args:
matrix: ndarray of size [m,n].
zero_tol: the tolerance to determine if two real number are identical or not.
Returns:
List[int]: A list of indices corresponding to the rows that are
lexicographical minima (the indices may be non-unique).
"""
if matrix.size == 0:
raise ValueError("The input matrix is empty")
if matrix.ndim > 2:
raise ValueError("The matrix must be 2 dimension")
if matrix.ndim == 1:
matrix = np.reshape(matrix, (1, matrix.size))
if matrix.shape[0] == 1:
return [0]
# Compare through a column by column manner. If the first element in the nth row is not the minimum, then we don't
# need to care about this row anymore because this row cannot be the lexico minimum.
[m, n] = matrix.shape
undetermined_rows = np.arange(m)
for col in range(n):
# find the minimum element in this column among all undetermined rows.
this_col = matrix[undetermined_rows, col]
argmin_in_this_col = np.where(abs(this_col - np.min(this_col)) < zero_tol)[0]
undetermined_rows = undetermined_rows[argmin_in_this_col]
if len(undetermined_rows) == 1:
break
return undetermined_rows.tolist()
[docs]
def is_symmetric(A: Matrix, zero_tol: float = tolerances['zero']) -> bool:
"""Check if a real matrix is symmetric or not."""
return bool(np.all((A - A.T) < zero_tol))
[docs]
def is_positive_definite(A: Matrix) -> bool:
r"""Check if a real matrix is positive definite (P.D.) through Cholesky decomposition.
Since :math:`x^\top Ax = (x^\top Ax)^\top = x^\top A^\top x, x^\top (A+A^\top)x = 2x^\top Ax`. Therefore, :math:`A` is P.D. if and only if :math:`A^\top+A` is P.D.
"""
try:
np.linalg.cholesky(A+A.T)
return True
except np.linalg.LinAlgError:
return False
[docs]
def is_positive_semidefinite(A: Matrix, hermitian: bool = False) -> bool:
r"""
Check if a real matrix is positive semi-definite (P.S.D.) through Cholesky decomposition.
Since :math:`x^\top Ax = (x^\top Ax)^\top = x^\top A^\top x, x^\top (A+A^\top)x = 2x^\top Ax`. Therefore, :math:`A` is P.S.D. if and only if :math:`A^\top+A` is P.S.D.
"""
if A.ndim != 2 or A.shape[0] != A.shape[1]:
return False
# Check if the matrix is symmetric
if not np.allclose(A, A.T):
return False
eigenvalues = np.linalg.eigvalsh(A)
return bool(np.all(eigenvalues >= 0))
[docs]
def is_p_matrix(A: Matrix) -> bool:
r"""Check if a matrix is a P-matrix.
A matrix is called P-matrix if its all principal minors are strictly positive. Refer to https://www.tuhh.de/ti3/paper/rump/Rump03b.pdf.
"""
if A.shape[0] != A.shape[1]:
raise ValueError("Matrix must be square")
n = A.shape[0]
# a positive definite matrix must be a P matrix
if is_positive_definite(A):
return True
else:
# iterate over all possible sizes of principal minors (1x1 to nxn)
for k in range(1, n+1):
# get all combinations of k rows/columns to form kxk submatrices
for rows_cols in combinations(range(n), k):
# extract the submatrix corresponding to the selected rows and columns
submatrix = A[np.ix_(rows_cols, rows_cols)]
if np.linalg.det(submatrix) <= 0:
return False
return True
[docs]
class Basis:
r"""
Define the basis of a Linear Complementarity Problem (LCP)
Any set :math:`B` which is a subset of :math:`\{1, ..., 2n\}` such that :math:`|B| = n` and :math:`\text{rank} A_{*,B} = n` is called a basis.
"""
def __init__(self, indices: List[int], n: Optional[int] = None):
# check if there are repeated elements in the indices list
if len(set(indices)) < len(indices):
raise ValueError("There are repeated elements in indices")
self._indices = indices.copy()
if n is None:
self._n = len(indices)
else:
self._n = n
if np.any(np.array(indices) < 0) or np.any(np.array(indices) > 2*self.n):
raise ValueError("Invalid value in indices")
@property
def indices(self):
return self._indices
@property
def n(self):
return self._n
def __eq__(self, other):
if isinstance(other, Basis):
return self.n == other.n and set(self.indices) == set(other.indices)
return False
def __len__(self):
return len(self.indices)
def __hash__(self):
# Combine the hash of the n attribute and the hash of the frozenset of indices
return hash((self.n, frozenset(self.indices)))
def __str__(self):
return f"Basis of length {self.n} with indices: {self.indices}"
def __repr__(self):
# return f"Basis of length {self.n} with indices: {self.indices}"
return str(self.indices)
[docs]
def is_complementary(self) -> bool:
"""
Determine if the current basis B is a complementary basis,
i.e., for all i in B, i+n and i-n are not in B
"""
for i in range(self.n):
if (self.indices[i] + self.n) % (2*self.n) in self.indices:
return False
return True
[docs]
def is_almost_complementary(self) -> bool:
"""
Determine if the current basis is almost complementary
"""
num_non_complementary = 0
for index in self.indices:
if (index + self.n) % (2*self.n) in self.indices:
num_non_complementary += 1
return True if num_non_complementary == 1 else False
[docs]
def is_feasible(self, M: Matrix, q: Vector) -> bool:
"""
Determine if the basis is feasible regrading a LCP (M, q)
"""
if M.ndim != 2 or M.shape[0] != M.shape[1]:
raise ValueError("M must be a sqaure matrix")
if q.ndim != 1:
raise ValueError("q must be a vector")
if M.shape[0] != q.size:
raise ValueError("Dimension mismatch between M and q")
A = np.hstack((np.eye(self.n), -M))
if np.all(np.linalg.solve(A[:, self.indices], q) >= -tolerances['zero']):
return True
else:
return False