import functools
from textwrap import indent
import numpy as np
from scipy.linalg import orth
from ddg.abc import LinearTransformable, NonExact
import ddg.math.projective as pmath
from ddg.math.euclidean import embed
import ddg.math.inner_product as ip
import ddg.math.linalg as la
[docs]class Subspace(LinearTransformable, NonExact):
"""Subspace of a projective space.
Represented by a number of vectors in homogeneous coordinates.
Parameters
----------
*points : array_like of list or numpy.ndarray of shape (n,)
Homogeneous coordinate vectors spanning the subspace. If the span is
1-dimensional, automatically casts to subspaces.Point
atol : float
rtol : float
Attributes
----------
points
matrix
dimension
ambient_dimension
codimension
Methods
-------
affine_points
affine_matrix
affine_point_and_directions
at_infinity
See Also
--------
Point
Notes
-----
Subspaces can be empty. In this case, points are stored as an empty matrix
with shape (n,0) if the ambient dimension is n-1. Empty spaces must be
initialized with 0-vectors.
"""
def __init__(self, *points, atol=None, rtol=None):
NonExact.__init__(self, atol=atol, rtol=rtol)
#Exception blocks
if len(np.shape(points[0])) != 1:
raise ValueError('Vectors spanning the subspace must be 1D '
'arrays given as separate arguments.')
if np.size(points[0]) == 0:
raise ValueError("Can't determine ambient dimension from empty "
"coordinate vectors.")
for p in points[1:]:
if np.shape(p) != np.shape(points[0]):
raise ValueError('Point vector length mismatch.')
point_matrix = np.array(points).T
# Reduce linearly dependent vectors.
# Also converts np.zeros(n) to an empty array of shape (n,0).
temp = la.col_basis(point_matrix, atol=self.atol, rtol=self.rtol)
if np.shape(temp) != np.shape(point_matrix):
point_matrix = temp
self._matrix = point_matrix
# Cast to Point. This would bypass Point.__init__ if it existed.
if np.shape(self._matrix)[1] == 1:
self.__class__ = Point
LinearTransformable.__init__(self, self.ambient_dimension+1)
@property
def points(self):
"""List of points spanning the subspace in homogeneous coordinates.
Returns
-------
list of numpy.ndarray of shape (ambient_dimension+1,)
"""
return la.matrix_to_arrays(self.matrix)
@property
def matrix(self):
"""Matrix whose columns span the subspace in homogeneous coordinates.
Returns
-------
numpy.ndarray of shape (ambient_dimension+1, dimension+1)
Raises
------
ValueError
If the subspace is at infinity.
"""
trafo = self.transformation
return np.dot(trafo, self._matrix)
[docs] def affine_points(self, affine_component=-1):
"""List of points spanning the affine part of the subspace.
Computes points in affine coordinates that, when homogenized according
to `affine_component`, form a basis for the subspace.
Parameters
----------
affine_component : int (default=-1)
Returns
-------
list of numpy.ndarray of shape (ambient_dimension,)
"""
return la.matrix_to_arrays(self.affine_matrix(affine_component))
[docs] def affine_matrix(self, affine_component=-1):
"""Matrix whose columns span the affine part of the subspace.
Computes points in affine coordinates that, when homogenized according
to `affine_component`, form a basis for the subspace. Returns them as
the columns of a matrix.
Parameters
----------
affine_component : int (default=-1)
Returns
-------
numpy.ndarray of shape (ambient_dimension, dimension+1)
Raises
------
ValueError
If the subspace is at infinity.
"""
if self.at_infinity(affine_component):
raise ValueError("Subspace at infinity.")
points = self.matrix
idx = affine_component
basemask = np.isclose(points[idx], 0,
atol=self.atol, rtol=self.rtol)
pointmask = np.where(basemask == False)[0]
directmask = np.where(basemask == True)[0]
points[:, pointmask] = points[:, pointmask]/ points[idx, pointmask]
point = points.T[pointmask[0]]
point = point[:, np.newaxis]
norms = np.linalg.norm(points[:, directmask], axis=0)
points[:, directmask] = points[:, directmask] / norms
points[:, directmask] = points[:, directmask] + point
return np.delete(points, idx, axis=0)
[docs] def affine_point_and_directions(self, affine_component=-1):
"""Return point in the subspace and directions.
Returns affine coordinate vectors p and d1, ..., dk such that
the original subspace is spanned by [p,1], [d1,0], ..., [dk,0]
where 1 or 0 are at position affine_component.
The point is the one that is closest to the zero vector in affine
coordinates and the directions are orthonormal.
Parameters
----------
affine_component : int (default=-1)
Returns
-------
numpy.ndarray, list of numpy.ndarray
See also
--------
subspace_from_affine_points_and_directions
"""
if self.at_infinity(affine_component):
raise ValueError("Subspace at infinity.")
points = self.matrix
idx = affine_component
basemask = np.isclose(points[idx], 0, atol=self.atol, rtol=self.rtol)
# dehomogenize the points that are not at infinity
pointmask = np.where(basemask == False)[0]
points[:, pointmask] = points[:, pointmask] / points[idx, pointmask]
# isolate the point to be returned
pidx = pointmask[0]
pointmask = pointmask[1:]
point = points.T[pidx]
point = point[:, np.newaxis]
# Directions as differences. The other directions are points at
# infinity and we can just delete the 0.
points[:, pointmask] -= point
if pidx != 0:
points[:, (0, pidx)] = points[:, (pidx, 0)]
temp = np.delete(points, idx, axis=0)
directions = temp[:, 1:]
if np.size(directions) != 0:
directions = orth(directions)
directions = list(directions.T)
point = temp[:, 0]
# project so the point is the closest one to the origin
point -= sum(np.dot(point, d)*d for d in directions)
return point, directions
[docs] def at_infinity(self, affine_component=-1):
"""Whether or not the subspace is at infinity.
Returns True if and only if p[affine_component] == 0 for all points p.
Parameters
----------
affine_component : int (default=-1)
Returns
-------
bool
"""
points = self.matrix
idx = affine_component
tmp = np.isclose(points[idx], 0, atol=self.atol, rtol=self.rtol)
return tmp.all()
@property
def dimension(self):
return np.shape(self._matrix)[1] - 1
@property
def ambient_dimension(self):
return np.shape(self._matrix)[0] - 1
@property
def codimension(self):
return self.ambient_dimension - self.dimension
[docs] def dualize(self, ambient_space=None):
"""Return dual subspace.
Parameters
----------
ambient_space : Subspace (default=None)
Other subspace containing this subspace that should be used as
ambient space when dualizing. None means whole space.
Returns
-------
Subspace
"""
# Probably not necessary, but we avoid an unnecessary meet
if ambient_space is None:
return subspace_from_columns(ip.get_col_complement(self.matrix))
if not (self <= ambient_space):
raise ValueError("Subspace must be contained in ambient space.")
return meet(self.dualize(), ambient_space)
[docs] def orthonormalize(self, affine_component=-1):
"""
Parameters
----------
affine_component : int (default=-1)
Returns
-------
Subspace
See also
--------
orthonormalize_subspace
"""
return orthonormalize_subspace(self, affine_component=affine_component)
[docs] def center(self, center, remove_index=None, insert_index=0,
affine_component=-1):
"""
Parameters
----------
affine_component : int (default=-1)
Returns
-------
Subspace
See also
--------
center_subspace
"""
return center_subspace(self, center, remove_index=remove_index,
insert_index=insert_index,
affine_component=affine_component)
[docs] def dehomogenize(self, affine_component=-1):
"""
Parameters
----------
affine_component : int (default=-1)
Returns
-------
Subspace
See also
--------
dehomogenize_subspace
"""
return dehomogenize_subspace(self, affine_component=affine_component)
def __contains__(self, p):
"""
Checks whether a homogeneous coordinate vector or Point object is
contained in the subspace.
"""
if isinstance(p, Point):
p = p.point
if len(p) != self.ambient_dimension + 1:
raise ValueError(f"Dimension mismatch: Neither a homogeneous "
f"coordinate vector nor Point object in "
f"{self.ambient_dimension}D space.")
eps = [np.dot(p, q) for q in self.dualize().points]
return all(np.isclose(e, 0., atol=self.atol, rtol=self.rtol)
for e in eps)
def __ge__(self, other):
"""
Checks if subspace contains `other`.
"""
if isinstance(other, Subspace):
if self.ambient_dimension != other.ambient_dimension:
raise ValueError(f"Ambient dimension mismatch: "
f"{self.ambient_dimension} and "
f"{other.ambient_dimension}.")
points = other.matrix
dual_matrix = self.dualize().matrix.T
eps = dual_matrix.dot(points)
return np.isclose(eps, np.zeros(eps.shape), atol=self.atol,
rtol=self.rtol).all()
else:
raise TypeError(repr(other) + ' is not of Subspace class.')
def __le__(self, other):
"""
Checks if subspace is contained in `other`.
"""
return other.__ge__(self)
def __eq__(self, other):
if not isinstance(other, Subspace):
return False
return (self <= other and self >= other)
def __repr__(self):
if self.dimension != -1:
points_param = ", ".join([repr(p) for p in self.points])
else:
points_param = repr(np.zeros(self.ambient_dimension + 1))
return f"Subspace({points_param})"
def __str__(self):
if self.dimension == -1:
return f"Empty space in {self.ambient_dimension}D projective space"
elif self.codimension == 0:
s = "Whole space"
elif self.dimension == 1:
s = "Line"
elif self.dimension == 2:
s = "Plane"
elif self.codimension == 1:
s = "Hyperplane"
else:
s = f"{self.dimension}D Subspace"
return (f"{s} in {self.ambient_dimension}D projective space\n"
f"Homogeneous coordinates (columns):\n"\
+ indent(str(self.matrix), ' '))
[docs]class Point(Subspace):
"""Subclass for points (0-dim. projective subspaces).
Parameters
----------
point : list or numpy.ndarray of shape (n,)
atol : float
rtol : float
Attributes
----------
point
affine_point
See also
--------
ddg.geometry.subspaces.Subspace
"""
def __init__(self, point, atol=None, rtol=None):
super().__init__(point, atol=atol, rtol=rtol)
@property
def point(self):
"""Homogeneous coordinate vector representing the point.
Returns
-------
numpy.ndarray (dimension+1,)
"""
return self.points[0]
[docs] def affine_point(self, affine_component=-1):
"""Affine coordinate vector representing the point.
Parameters
----------
affine_component : int (default=-1)
Returns
-------
numpy.ndarray of shape (dimension,)
Raises
------
ValueError
If the point is at infinity.
"""
return self.affine_points(affine_component)[0]
def __repr__(self):
return f"Point({repr(self.point)})"
def __str__(self):
return (f"Point in {self.ambient_dimension}D projective space\n"
f"Homogeneous coordinates: {self.point}")
[docs]def subspace_from_rows(matrix, **kwargs):
"""Create subspace using the rows of a matrix as points.
Parameters
----------
matrix : numpy.ndarray of shape (k+1,n+1)
**kwargs : dict
Other arguments to be passed to the Subspace.
Returns
-------
ddg.geometry.subspaces.Subspace
Notes
-----
An empty subspace can be created by giving a 0-matrix or a matrix with
shape (0,n+1).
"""
if np.size(matrix) == 0:
matrix = [np.zeros(np.shape(matrix)[1])]
return Subspace(*matrix, **kwargs)
[docs]def subspace_from_columns(matrix, **kwargs):
"""Create subspace using the columns of a matrix as points.
Parameters
----------
matrix : numpy.ndarray of shape (n+1,k+1)
**kwargs : dict
Other arguments to be passed to the Subspace.
Returns
-------
ddg.geometry.subspaces.Subspace
Notes
-----
An empty subspace can be created by giving a 0-matrix or a matrix with
shape (n,0).
"""
return subspace_from_rows(matrix.T, **kwargs)
[docs]def hyperplane_from_normal(normal, level=None, point=None, **kwargs):
"""
An affine hyperplane defined by either
<x, normal> == level,
or the normal and a point (given in affine coordinates) on the plane.
"""
if (level is None) == (point is None):
raise ValueError("Hyperplane must be specified by a normal and either "
"the level or one point.")
elif level is not None:
h = Point(np.concatenate([normal, [-level]]), **kwargs).dualize()
else:
level = np.dot(normal, point)
h = hyperplane_from_normal(normal, level=level, **kwargs)
if 'atol' in kwargs:
h.atol = kwargs['atol']
if 'rtol' in kwargs:
h.rtol = kwargs['rtol']
return h
[docs]def normals(subspace, affine_component=-1):
"""Return all normal directions of a subspace.
If p1, p2 are two points in the subspace in affine coordinates and n is a
normal direction, then
<p1 - p2, n> = 0.
Parameters
----------
subspace : Subspace
affine_component : int (default=-1)
Returns
-------
list of numpy.ndarray of shape (subspace.ambient_dimension,)
Normals in affine coordinates.
"""
i = affine_component
if subspace.dimension < 1:
return whole_space(subspace.ambient_dimension)
dirs = subspace.affine_point_and_directions(i)[1]
dirs = np.column_stack(dirs)
normals = ip.get_col_complement(dirs)
return list(normals.T)
[docs]def normal_with_level(hyperplane, affine_component=-1):
if hyperplane.codimension != 1:
raise TypeError("Only hyperplanes have normals.")
if hyperplane.at_infinity(affine_component):
raise ValueError("Hyperplane is at infinity.")
n = np.ma.array(hyperplane.dualize().point, mask=False)
ac = affine_component
n.mask[ac] = True
norm = np.linalg.norm(n.compressed())
n.mask[ac] = False
n = (n/norm).compressed()
if not ac in (-1, hyperplane.ambient_dimension +1):
n[ac], n[-1] = n[-1], n[ac]
return n
[docs]def normal(hyperplane, affine_component=-1):
return normal_with_level(hyperplane, affine_component)[:-1]
[docs]def level(hyperplane, affine_component=-1):
return -normal_with_level(hyperplane, affine_component)[-1]
[docs]def subspace_from_affine_points(*points, affine_component=-1, **kwargs):
"""Create a subspace from affine points.
Parameters
----------
*points : iterable of numpy.ndarray or list
Affine coordinate vectors
affine_component : int (optional, default=-1)
Where to insert the affine_component for the subspace. The default -1
will append to the end.
**kwargs : dict
Other keyword arguments to be passed to the subspace.
Returns
-------
ddg.geometry.subspaces.Subspace
Subspace spanned by the given points
"""
points = [pmath.homogenize(p, affine_component) for p in points]
return Subspace(*points, **kwargs)
[docs]def subspace_from_affine_rows(matrix, affine_component=-1, **kwargs):
"""Create a subspace from the rows of a matrix, interpreted as affine
coordinates.
Parameters
----------
matrix : numpy.ndarray
Array whichs columns describe the affine points
affine_component : int (optional, default=-1)
Where to insert the affine_component for the subspace. The default -1
will append to the end.
**kwargs : dict
Other keyword arguments to be passed to the subspace.
Returns
-------
ddg.geometry.subspaces.Subspace
Subspace spanned by the given points
"""
return subspace_from_affine_points(*matrix,
affine_component=affine_component,
**kwargs)
[docs]def subspace_from_affine_columns(matrix, affine_component=-1, **kwargs):
"""Create a subspace from the columns of a matrix, interpreted as affine
coordinates.
Parameters
----------
matrix : numpy.ndarray
Array whichs columns describe the affine points
affine_component : int (optional, default=-1)
Where to insert the affine_component for the subspace. The default -1
will append to the end.
**kwargs : dict
Other keyword arguments to be passed to the subspace.
Returns
-------
ddg.geometry.subspaces.Subspace
Subspace spanned by the given points
"""
return subspace_from_affine_rows(matrix.T,
affine_component=affine_component,
**kwargs)
[docs]def subspace_from_affine_points_and_directions(points, directions,
affine_component=-1, **kwargs):
"""Create a subspace from affine points and directions.
Parameters
----------
points : (list of numpy.array) or numpy.array
List of points or a single point
directions : (list of numpy.array) or numpy.array
List of directions or a single direction
affine_component : int (optional, default=-1)
Where to insert the affine_component for the subspace. The default -1
will append to the end.
**kwargs : dict
Other keyword arguments to be passed to the subspace.
Returns
-------
ddg.geometry.subspaces.Subspace
Subspace spanned by the given points and directions
Raises
------
ValueError
If both points and directions are empty lists.
"""
if len(np.shape(points)) == 1 and len(points) != 0:
points = [points]
if len(np.shape(directions)) == 1 and len(directions) != 0:
directions = [directions]
num_pts = len(points)
pts_and_dirs = points + directions
if len(pts_and_dirs) == 0:
raise ValueError("Points and directions can not both be empty.")
for i, p in enumerate(pts_and_dirs[:num_pts]):
pts_and_dirs[i] = pmath.homogenize(
p,
affine_component=affine_component
)
for i, d in enumerate(pts_and_dirs[num_pts:]):
pts_and_dirs[num_pts + i] = embed(
d,
affine_component=affine_component
)
return Subspace(*pts_and_dirs, **kwargs)
[docs]def orthonormalize_subspace(subspace, affine_component=-1):
"""Regularize the basis of a subspace.
Takes a subspace and returns another subspace with a regularized basis
based on ``subspace.affine_point_and_directions()``. Useful in preparation
for visualization.
Parameters
----------
subspace : Subspace
affine_component : int (default=-1)
Returns
-------
Subspace
"""
return subspace_from_affine_points_and_directions(
*subspace.affine_point_and_directions(affine_component),
affine_component=affine_component,
atol=subspace.atol,
rtol=subspace.rtol
)
[docs]def center_subspace(subspace, center, remove_index=None, insert_index=0,
affine_component=-1):
"""Insert a specific point into the basis of a subspace.
This is especially useful in preparation for visualization: When applying
this function and then using the "convex" parametrization in
`to_smooth_net`, the parametrization will be centered around the new
`center`, if it is a point not at infinity.
Parameters
----------
subspace : Subspace
center : Point or numpy.ndarray of shape (k,) or (k+1,)
Where k is the ambient dimension of the subspace. Point must lie in
subspace. Can be given as Point object or as array containing affine or
homogeneous coordinates.
remove_index : int (default=None)
Index of vector that should be removed in favor of center. If None, we
search for one.
insert_index : int (default=0)
Index where center should be inserted.
affine_component : int (default=-1)
Returns
-------
Subspace
The same subspace with a basis that has `center` at position
insert_index.
Raises
------
ValueError
* If center does not lie in the subspace.
* If the vector at `remove_index` could not be replaced with `center`
while still having a basis.
"""
if isinstance(center, Point):
center = center.point
elif np.size(center) == subspace.ambient_dimension:
center = pmath.homogenize(center, affine_component=affine_component)
if center not in subspace:
raise ValueError("Center must lie in subspace.")
# Let v1,...,vk be a basis. If w = a1 * v1 + ... + ak * vk and ai != 0,
# then v1,...,w,...,vk is again a basis, where vi was exchanged with w.
coords = la.coordinates(center, subspace.matrix, atol=subspace.atol,
rtol=subspace.rtol)
if remove_index is not None:
idx = remove_index
if np.isclose(coords[idx], 0, atol=subspace.atol, rtol=subspace.rtol):
raise ValueError(
f"Could not replace vector {subspace.points[remove_index]} "
f"with center {center} and still have a basis."
)
else:
for idx, coord in enumerate(coords):
if not np.isclose(coord, 0, atol=subspace.atol,
rtol=subspace.rtol):
break
mask = np.full(subspace.dimension + 1, True)
mask[idx] = False
# Convert to non-negative index, because otherwise insert will act
# strangely.
if -(subspace.dimension + 1) <= insert_index <= -1:
insert_index += (subspace.dimension + 1)
return subspace_from_columns(
np.insert(subspace.matrix[:, mask], insert_index, center, axis=1),
atol=subspace.atol,
rtol=subspace.rtol
)
[docs]def dehomogenize_subspace(subspace, affine_component=-1):
"""Normalize affine_component entries to 1 if possible.
Useful for visualization.
Parameters
----------
Subspace
affine_component : int (default=-1)
Returns
-------
Subspace
"""
entries = subspace.matrix[affine_component]
diag = []
for e in entries:
if np.isclose(e, 0, atol=subspace.atol, rtol=subspace.rtol):
diag.append(1)
else:
diag.append(1 / e)
return subspace_from_columns(
subspace.matrix @ np.diag(diag),
atol=subspace.atol,
rtol=subspace.rtol
)
[docs]def whole_space(dimension, **kwargs):
"""Create whole space of a certain dimension as a trivial subspace.
Parameters
----------
Dimension : int
Dimension of the space.
**kwargs
Keyword arguments other than points to be passed to Subspace.
Returns
-------
ddg.geometry.subspaces.Subspace
See also
--------
Subspace
"""
return subspace_from_columns(np.eye(dimension+1), **kwargs)
@NonExact.nonexact_function
def _join(subspace1, subspace2, atol=None, rtol=None):
"""Construct the join of two subspaces.
Parameters
----------
subspace1 : ddg.geometry.subspaces.Subspace
Subspace to join with subspace2
subspace2 : ddg.geometry.subspaces.Subspace
Subspace to join with subspace1
atol : float (default=None)
rtol : float (default=None)
Returns
-------
ddg.geometry.subspaces.Subspace
Raises
------
ValueError
If ambient_dimension of the subspaces does not match.
See Also
--------
meet, intersect
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See ddg.abc.NonExact for details.
"""
if subspace1.ambient_dimension != subspace2.ambient_dimension:
raise ValueError(f"Ambient dimension mismatch: "
f"{subspace1.ambient_dimension} and "
f"{subspace2.ambient_dimension}.")
return Subspace(
*(subspace1.points + subspace2.points),
atol=atol,
rtol=rtol,
)
[docs]def join(*subspaces):
"""Construct the join of linear subspaces.
Parameters
----------
*subspaces : ddg.geometry.subspaces.Subspace
Returns
-------
result : ddg.geometry.subspaces.Subspace
Join of the subspaces
Notes
-----
For the tolerances, we use the maximum of the tolerances of the subspaces.
"""
atol = max([s.atol for s in subspaces])
rtol = max([s.rtol for s in subspaces])
_join_with_tols = functools.partial(_join, atol=atol, rtol=rtol)
return functools.reduce(_join_with_tols, subspaces)
@NonExact.nonexact_function
def _meet(subspace1, subspace2, atol=None, rtol=None):
"""Construct the meet of two linear subspaces.
Parameters
----------
subspace1 : ddg.geometry.subspaces.Subspace
Subspace to meet with subspace2
subspace2 : ddg.geometry.subspaces.Subspace
Subspace to meet with subspace1
atol : float (default=None)
rtol : float (default=None)
Returns
-------
ddg.geometry.subspaces.Subspace
Meet of the subspaces
Raises
------
ValueError
If ambient_dimension of the subspaces does not match.
See Also
--------
join, intersect
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See ddg.abc.NonExact for details.
"""
if subspace1.ambient_dimension != subspace2.ambient_dimension:
raise ValueError(f"Ambient dimension mismatch: "
f"{subspace1.ambient_dimension} and "
f"{subspace2.ambient_dimension}.")
points1 = subspace1.dualize().matrix
points2 = subspace2.dualize().matrix
points = np.append(points1, points2, axis=1)
points = la.col_basis(points, atol=atol, rtol=rtol)
result = subspace_from_columns(points).dualize()
return result
[docs]def meet(*subspaces):
"""Construct the meet of projective subspaces.
Parameters
----------
*subspaces : ddg.geometry.subspaces.Subspace
Returns
-------
result : ddg.geometry.subspaces.Subspace
Meet of the subspaces
Notes
-----
For the tolerances, we use the maximum of the tolerances of the subspaces.
"""
atol = max([s.atol for s in subspaces])
rtol = max([s.rtol for s in subspaces])
_meet_with_tols = functools.partial(_meet, atol=atol, rtol=rtol)
return functools.reduce(_meet_with_tols, subspaces)
[docs]def intersect(*subspaces):
"""Construct intersection of projective subspaces.
This is an alias for meet.
See Also
--------
meet, join
"""
return meet(*subspaces)
[docs]def least_square_subspace(points, k, affine_component=-1):
"""
Computes the k-dimensional subspace closest to N projective
points. See below for more details.
Parameters
----------
points : list of ddg.Point or numpy.ndarray, each of shape (n+1,)
The given N projective points.
k : int
Required dimension of the subspace.
affine_component : int (default=-1)
Index of affine component normalized to one.
Returns
-------
subspace : ddg.Subspace
k-dimensional projective subspace.
Raises
------
ValueError
if no points are given.
TypeError
if points are not given as a list of ddg.Points
nor as a list of numpy.ndarray.
ValueError
if the given points (of type numpy.ndarray) are not
one-dimensional numpy.ndarray.
ValueError
if k is strictly greater than the ambient projective
dimension n.
ValueError
if the points do not all lie in the n-dimensional copy of
affine space defined by x_j != 0 where j is the coordinate
given by affine_component.
Notes
-----
The resulting subspace is 'best-fitting' or 'closest' k-dimensional
subspace to the N points in the sense that it minimizes the sum of
the squares of the Euclidean distances of the points to it.
For an appropriate k << N, the subspace is computed in the following
steps:
- First, dehomogenize the N points and define a (N x n)-matrix
'A' whoes rows are their chosen affine coordinates.
- Second, compute the matrix 'B' with rank k+1 that is closest
to A in the Frobenius norm. This is done by first decomposing A
using SVD into A = U * np.diag(s) * Vh where U and Vh are unitary
matrices of shape (N x N) and (n x n) respectively and s is 1D
array of singular values of A. Then defining B as the sum of the
(N x n)-matrices, each of rank 1, given by product '(u * s) * vh'
truncated after the first k+1 terms.
- Finally, return the projectivization of the (k+1)-dimensional
subspace defined by the matrix B.
"""
if isinstance(points, list):
if points:
__points = []
for point in points:
if isinstance(point, Point):
__points.append(point)
elif isinstance(point, np.ndarray):
if point.ndim == 1:
__points.append(Point(point))
else:
raise ValueError(f"Any point given as numpy.ndarray must "
f"be 1D, not {point.ndim}D numpy.ndarray.")
else:
raise TypeError(f"All given points must be of type Point or "
f"numpy.ndarray, not {type(point).__name__}.")
else:
raise ValueError("No points are given to return least-square subspace.")
points_span = Subspace(*[p.point for p in __points])
n = points_span.ambient_dimension
r = points_span.dimension
if k > n:
raise ValueError("Can not return subspace with dimension greater than"
" ambient dimension.")
elif k == n:
subspace = whole_space(k)
elif k > r:
basis_matrix = la.nullspace(points_span.matrix.T)
subspace = subspace_from_columns(np.column_stack((points_span.matrix,
basis_matrix[:, 0:(k - r)])))
elif k == r:
subspace = points_span
else:
if all(q.point[affine_component] != 0 for q in __points):
affine_coord = [pmath.dehomogenize(q.point, affine_component) for q in __points]
else:
raise ValueError(f"The given points do not lie in a common copy of"
f" {n}-dimensional affine space.")
A = np.array(affine_coord)
U, S, Vh = np.linalg.svd(A)
u = U[:, :k]
s = np.diag(S)[:k, :k]
vh = Vh[:, :k].T
A_kth_truncate = (u @ s) @ vh
subspace = Subspace(*[pmath.homogenize(np.ravel(A_kth_truncate[i, :]))
for i in range(0, A_kth_truncate.shape[0])])
return subspace
else:
type_error_msg = f"{points.ndim}D numpy.ndarray" if isinstance(points, np.ndarray) \
else f"{type(points).__name__}"
raise TypeError(f"Input points must given as a list of Points or a list of "
f"1D numpy.ndarray, not as {type_error_msg}.")