Find maximal angles between polyhedral convex cones#
Warning
This module is considered internal and its contents are subject to
change at any time without (deprecation) warning. The stable
interface is
sage.geometry.cone.ConvexRationalPolyhedralCone.max_angle()
.
Finding the maximal (or equivalently, the minimal) angle between two polyhedral convex cones is a hard nonconvex optimization problem. The problem for a single cone was introduced in [IS2005], and was later extended in [SS2016] to two cones as a generalization of the principal angle between two vector subspaces.
Seeger and Sossa proposed an algorithm in [SS2016] to find maximal angles, and [Or2020] elaborates on that algorithm. It is this latest improvement that is implemented (more or less) by this module. The fact that perturbations of pointed cones may not change the answer too much [Or2024] is taken into consideration to avoid pathological cases.
This module is internal to SageMath; the interface presented to users
consists of a public method,
sage.geometry.cone.ConvexRationalPolyhedralCone.max_angle()
for
polyhedral convex cones. Even though all of the functions in this
module are internal, some are more internal than others. There are a
few functions that are used only in doctests, and not by any code that
an end-user would run. Breaking somewhat with tradition, only those
methods have been prefixed with an underscore.
- sage.geometry.cone_critical_angles.check_gevp_feasibility(cos_theta, xi, eta, G_I, G_I_c_T, H_J, H_J_c_T, epsilon)[source]#
Determine if a solution to the generalized eigenvalue problem in Theorem 3 [Or2020] is feasible.
Implementation detail: we take four matrices that we are capable of computing as parameters instead, because we will be called in a nested loop “for all \(I\)… and for all \(J\)…” The data corresponding to \(I\) should be computed only once, which means that we can’t do it here – it needs to be done outside of the \(J\) loop. For symmetry (and to avoid relying on too many cross-function implementation details), we also insist that the \(J\) data be passed in.
INPUT:
cos_theta
– an eigenvalue corresponding to \(( \xi, \eta )\)xi
– first component of the \(( \xi, \eta )\) eigenvectoreta
– second component of the \(( \xi, \eta )\) eigenvectorG_I
– the submatrix of \(G\) with columns indexed by \(I\)G_I_c_T
– a matrix whose rows are the non-\(I\) columns of \(G\)H_J
– the submatrix of \(H\) with columns indexed by \(J\)H_J_c_T
– a matrix whose rows are the non-\(J\) columns of \(H\)epsilon
– the tolerance to use when making comparisons
OUTPUT:
A triple containing (in order),
a boolean,
a vector in the cone \(P\) (of the same length as
xi
), anda vector in the cone \(Q\) (of the same length as
eta
).
If \(( \xi, \eta )\) is feasible, we return
(True, u, v)
where \(u\) and \(v\) are the vectors in \(P\) and \(Q\) respectively that form the the angle \(\theta\).If \(( \xi, \eta )\) is not feasible, then we return
(False, 0, 0)
where0
should be interpreted to mean the zero vector in the appropriate space.EXAMPLES:
If \(\xi\) has any components less than “zero,” it isn’t feasible:
sage: from sage.geometry.cone_critical_angles import( ....: check_gevp_feasibility) sage: xi = vector(QQ, [-1,1]) sage: eta = vector(QQ, [1,1,1]) sage: check_gevp_feasibility(0,xi,eta,None,None,None,None,0) (False, (0, 0), (0, 0, 0))
>>> from sage.all import * >>> from sage.geometry.cone_critical_angles import( ... check_gevp_feasibility) >>> xi = vector(QQ, [-Integer(1),Integer(1)]) >>> eta = vector(QQ, [Integer(1),Integer(1),Integer(1)]) >>> check_gevp_feasibility(Integer(0),xi,eta,None,None,None,None,Integer(0)) (False, (0, 0), (0, 0, 0))
If \(\eta\) has any components less than “zero,” it isn’t feasible:
sage: from sage.geometry.cone_critical_angles import( ....: check_gevp_feasibility) sage: xi = vector(QQ, [2]) sage: eta = vector(QQ, [1,-4,4,5]) sage: check_gevp_feasibility(0,xi,eta,None,None,None,None,0) (False, (0), (0, 0, 0, 0))
>>> from sage.all import * >>> from sage.geometry.cone_critical_angles import( ... check_gevp_feasibility) >>> xi = vector(QQ, [Integer(2)]) >>> eta = vector(QQ, [Integer(1),-Integer(4),Integer(4),Integer(5)]) >>> check_gevp_feasibility(Integer(0),xi,eta,None,None,None,None,Integer(0)) (False, (0), (0, 0, 0, 0))
If \(\xi\) and \(\eta\) are equal and if \(G_{I}\) and \(H_{J}\) are not, then the copy of \(\eta\) that’s been scaled by the norm of \(G_{I}\xi\) generally won’t satisfy its norm-equality constraint:
sage: from sage.geometry.cone_critical_angles import( ....: check_gevp_feasibility) sage: xi = vector(QQ, [1,1]) sage: eta = xi sage: G_I = matrix.identity(QQ,2) sage: H_J = 2*G_I sage: check_gevp_feasibility(0,xi,eta,G_I,None,H_J,None,0) (False, (0, 0), (0, 0))
>>> from sage.all import * >>> from sage.geometry.cone_critical_angles import( ... check_gevp_feasibility) >>> xi = vector(QQ, [Integer(1),Integer(1)]) >>> eta = xi >>> G_I = matrix.identity(QQ,Integer(2)) >>> H_J = Integer(2)*G_I >>> check_gevp_feasibility(Integer(0),xi,eta,G_I,None,H_J,None,Integer(0)) (False, (0, 0), (0, 0))
When \(\cos\theta\) is zero, the inequality (42) in Theorem 7.3 [SS2016] is just an inner product with \(v\) which we can make positive by ensuring that all of the entries of \(H_{J}\) are positive. So, if any of the rows of
G_I_c_T
contain a negative entry, (42) will fail:sage: from sage.geometry.cone_critical_angles import( ....: check_gevp_feasibility) sage: xi = vector(QQ, [1/2,1/2,1/2,1/2]) sage: eta = xi sage: G_I = matrix.identity(QQ,4) sage: G_I_c_T = matrix(QQ, [[0,-1,0,0]]) sage: H_J = G_I sage: check_gevp_feasibility(0,xi,eta,G_I,G_I_c_T,H_J,None,0) (False, (0, 0, 0, 0), (0, 0, 0, 0))
>>> from sage.all import * >>> from sage.geometry.cone_critical_angles import( ... check_gevp_feasibility) >>> xi = vector(QQ, [Integer(1)/Integer(2),Integer(1)/Integer(2),Integer(1)/Integer(2),Integer(1)/Integer(2)]) >>> eta = xi >>> G_I = matrix.identity(QQ,Integer(4)) >>> G_I_c_T = matrix(QQ, [[Integer(0),-Integer(1),Integer(0),Integer(0)]]) >>> H_J = G_I >>> check_gevp_feasibility(Integer(0),xi,eta,G_I,G_I_c_T,H_J,None,Integer(0)) (False, (0, 0, 0, 0), (0, 0, 0, 0))
Likewise we can make (43) fail in exactly the same way:
sage: from sage.geometry.cone_critical_angles import( ....: check_gevp_feasibility) sage: xi = vector(QQ, [1/2,1/2,1/2,1/2]) sage: eta = xi sage: G_I = matrix.identity(QQ,4) sage: G_I_c_T = matrix(QQ, [[0,1,0,0]]) sage: H_J = G_I sage: H_J_c_T = matrix(QQ, [[0,-1,0,0]]) sage: check_gevp_feasibility(0,xi,eta,G_I,G_I_c_T,H_J,H_J_c_T,0) (False, (0, 0, 0, 0), (0, 0, 0, 0))
>>> from sage.all import * >>> from sage.geometry.cone_critical_angles import( ... check_gevp_feasibility) >>> xi = vector(QQ, [Integer(1)/Integer(2),Integer(1)/Integer(2),Integer(1)/Integer(2),Integer(1)/Integer(2)]) >>> eta = xi >>> G_I = matrix.identity(QQ,Integer(4)) >>> G_I_c_T = matrix(QQ, [[Integer(0),Integer(1),Integer(0),Integer(0)]]) >>> H_J = G_I >>> H_J_c_T = matrix(QQ, [[Integer(0),-Integer(1),Integer(0),Integer(0)]]) >>> check_gevp_feasibility(Integer(0),xi,eta,G_I,G_I_c_T,H_J,H_J_c_T,Integer(0)) (False, (0, 0, 0, 0), (0, 0, 0, 0))
Finally, if we ensure that everything works, we get back a feasible result along with the vectors (scaled \(\xi\) and \(\eta\)) that worked:
sage: from sage.geometry.cone_critical_angles import( ....: check_gevp_feasibility) sage: xi = vector(QQ, [1/2,1/2,1/2,1/2]) sage: eta = xi sage: G_I = matrix.identity(QQ,4) sage: G_I_c_T = matrix(QQ, [[0,1,0,0]]) sage: H_J = G_I sage: H_J_c_T = matrix(QQ, [[0,1,0,0]]) sage: check_gevp_feasibility(0,xi,eta,G_I,G_I_c_T,H_J,H_J_c_T,0) (True, (1/2, 1/2, 1/2, 1/2), (1/2, 1/2, 1/2, 1/2))
>>> from sage.all import * >>> from sage.geometry.cone_critical_angles import( ... check_gevp_feasibility) >>> xi = vector(QQ, [Integer(1)/Integer(2),Integer(1)/Integer(2),Integer(1)/Integer(2),Integer(1)/Integer(2)]) >>> eta = xi >>> G_I = matrix.identity(QQ,Integer(4)) >>> G_I_c_T = matrix(QQ, [[Integer(0),Integer(1),Integer(0),Integer(0)]]) >>> H_J = G_I >>> H_J_c_T = matrix(QQ, [[Integer(0),Integer(1),Integer(0),Integer(0)]]) >>> check_gevp_feasibility(Integer(0),xi,eta,G_I,G_I_c_T,H_J,H_J_c_T,Integer(0)) (True, (1/2, 1/2, 1/2, 1/2), (1/2, 1/2, 1/2, 1/2))
- sage.geometry.cone_critical_angles.compute_gevp_M(gs, hs)[source]#
Compute the matrix \(M\) whose \((i,j)\)-th entry is the inner product of
gs[i]
andhs[j]
.This is the “generalized Gram matrix” appearing in Proposition 6 in [Or2020]. For efficiency, we also return the minimal pair, whose inner product is minimal among the entries of \(M\). This allows our consumer to bail out immediately (knowing the optimal pair!) if it turns out that the maximal angle is acute; i.e. if the smallest entry of \(M\) is nonnegative.
INPUT:
gs
– a linearly independent list of unit-norm generators for the cone \(P\)hs
– a linearly independent list of unit-norm generators for the cone \(Q\)
OUTPUT:
A tuple containing four elements, in order:
The matrix \(M\) described in Proposition 6
The minimal entry in the matrix \(M\)
A vector in
gs
that achieves that minimal inner product along with the next element of the tupleA vector in
hs
that achieves the minimal inner product along with the previous element in the tuple
EXAMPLES:
sage: from sage.geometry.cone_critical_angles import compute_gevp_M sage: P = Cone([ (1,2,0), (3,4,0) ]) sage: Q = Cone([ (-1,4,1), (5,-2,-1), (-1,-1,5) ]) sage: gs = [g.change_ring(QQ) for g in P] sage: hs = [h.change_ring(QQ) for h in Q] sage: M = compute_gevp_M(gs, hs)[0] sage: all( M[i][j] == gs[i].inner_product(hs[j]) ....: for i in range(P.nrays()) ....: for j in range(Q.nrays()) ) True
>>> from sage.all import * >>> from sage.geometry.cone_critical_angles import compute_gevp_M >>> P = Cone([ (Integer(1),Integer(2),Integer(0)), (Integer(3),Integer(4),Integer(0)) ]) >>> Q = Cone([ (-Integer(1),Integer(4),Integer(1)), (Integer(5),-Integer(2),-Integer(1)), (-Integer(1),-Integer(1),Integer(5)) ]) >>> gs = [g.change_ring(QQ) for g in P] >>> hs = [h.change_ring(QQ) for h in Q] >>> M = compute_gevp_M(gs, hs)[Integer(0)] >>> all( M[i][j] == gs[i].inner_product(hs[j]) ... for i in range(P.nrays()) ... for j in range(Q.nrays()) ) True
- sage.geometry.cone_critical_angles.gevp_licis(G)[source]#
Return all nonempty subsets of indices for the columns of
G
that correspond to linearly independent sets (of columns ofG
).Mnemonic: linearly independent column-index subsets (LICIS).
The returned lists are all sorted in the same (the natural) order; and are returned as lists so that they may be used to index into the rows/columns of matrices.
INPUT:
G
– the matrix whose linearly independent column index sets we want
OUTPUT:
A generator that returns sorted lists of natural numbers. Each generated list
I
is a set of indices corresponding to columns ofG
that, when considered as a set, is linearly independent.EXAMPLES:
The linearly independent subsets of the matrix corresponding to a line (with two generators pointing in opposite directions) are the one-element subsets, since the only two-element subset isn’t linearly independent:
sage: from sage.geometry.cone_critical_angles import gevp_licis sage: K = Cone([(1,0),(-1,0)]) sage: G = matrix.column(K.rays()) sage: list(gevp_licis(G)) [[0], [1]]
>>> from sage.all import * >>> from sage.geometry.cone_critical_angles import gevp_licis >>> K = Cone([(Integer(1),Integer(0)),(-Integer(1),Integer(0))]) >>> G = matrix.column(K.rays()) >>> list(gevp_licis(G)) [[0], [1]]
The matrix for the trivial cone has no linearly independent subsets, since we require them to be nonempty:
sage: from sage.geometry.cone_critical_angles import gevp_licis sage: trivial_cone = cones.trivial(0) sage: trivial_cone.is_trivial() True sage: list(gevp_licis(matrix.column(trivial_cone.rays()))) []
>>> from sage.all import * >>> from sage.geometry.cone_critical_angles import gevp_licis >>> trivial_cone = cones.trivial(Integer(0)) >>> trivial_cone.is_trivial() True >>> list(gevp_licis(matrix.column(trivial_cone.rays()))) []
All rays in the nonnegative orthant of \(R^{n}\) are linearly independent, so we should get back \(2^{n} - 1\) subsets after accounting for the absence of the empty set:
sage: from sage.geometry.cone_critical_angles import gevp_licis sage: K = cones.nonnegative_orthant(3) sage: G = matrix.column(K.rays()) sage: len(list(gevp_licis(G))) == 2^(K.nrays()) - 1 True
>>> from sage.all import * >>> from sage.geometry.cone_critical_angles import gevp_licis >>> K = cones.nonnegative_orthant(Integer(3)) >>> G = matrix.column(K.rays()) >>> len(list(gevp_licis(G))) == Integer(2)**(K.nrays()) - Integer(1) True
- sage.geometry.cone_critical_angles.max_angle(P, Q, exact, epsilon)[source]#
Find the maximal angle between the cones \(P\) and \(Q\).
This implements
sage.geometry.cone.ConvexRationalPolyhedralCone.max_angle()
, which should be fully documented.EXAMPLES:
For the sake of the user interface, the argument validation for this function is performed in the associated cone method; we can therefore crash it by feeding it invalid input like an inadmissible cone:
sage: from sage.geometry.cone_critical_angles import max_angle sage: K = cones.trivial(3) sage: max_angle(K,K,True,0) Traceback (most recent call last): ... IndexError: list index out of range
>>> from sage.all import * >>> from sage.geometry.cone_critical_angles import max_angle >>> K = cones.trivial(Integer(3)) >>> max_angle(K,K,True,Integer(0)) Traceback (most recent call last): ... IndexError: list index out of range
- sage.geometry.cone_critical_angles.solve_gevp_nonzero(GG, HH, M, I, J)[source]#
Solve the generalized eigenvalue problem in Theorem 3 [Or2020] for a nonzero eigenvalue using Propositions 3 and 5 [Or2020].
INPUT:
GG
– the matrix whose \((i,j)\)-th entry is the inner product of \(g_{i}\) and \(g_{j}\), which are in turn the \(i\)-th and \(j\)-th columns of the matrix \(G\) in Theorem 3 [Or2020]HH
– the matrix whose \((i,j)\)-th entry is the inner product of \(h_{i}\) and \(h_{j}\), which are in turn the \(i\)-th and \(j\)-th columns of the matrix \(H\) in Theorem 3 [Or2020]M
– the matrix whose \((i,j)\)-th entry is the inner product of \(g_{i}\) and \(h_{j}\) as in Proposition 6 in [Or2020]I
– a linearly independent column-index set for the matrix \(G\) that appears in Theorem 3 [Or2020]J
– a linearly independent column-index set for the matrix \(H\) that appears in Theorem 3 [Or2020]
OUTPUT:
A generator of
(eigenvalue, xi, eta, multiplicity)
quartets whereeigenvalue
is a real eigenvalue of the systemxi
is the first (lengthlen(I)
) component of an eigenvector associated witheigenvalue
eta
is the second (lengthlen(J)
) component of an eigenvector associated witheigenvalue
multiplicity
is the dimension of the eigenspace associated witheigenvalue
Note that we do not return a basis for each eigenspace along with its eigenvalue. For the application we have in mind, an eigenspace of dimension greater than one (so,
multiplicity > 1
) is an error. As such, our return value is optimized for convenience in the non-error case, where there is only one eigenvector (spanning a one-dimensional eigenspace) associated with each eigenvalue.ALGORITHM:
According to Proposition 5 [Or2020], the solutions corresponding to non-zero eigenvalues can be found by solving a smaller eigenvalue problem in only the variable \(\xi\). So, we do that, and then solve for \(\eta\) in terms of \(\xi\) as described in the proposition.
EXAMPLES:
When the zero solutions are included, this function returns the same solutions as the naive method on the Schur cone in three dimensions:
sage: from itertools import chain sage: from sage.geometry.cone_critical_angles import ( ....: _normalize_gevp_solution, ....: _solve_gevp_naive, ....: gevp_licis, ....: solve_gevp_nonzero, ....: solve_gevp_zero) sage: K = cones.schur(3) sage: gs = [g.change_ring(AA).normalized() for g in K] sage: G = matrix.column(gs) sage: GG = G.transpose() * G sage: G_index_sets = list(gevp_licis(G)) sage: all( ....: set( ....: _normalize_gevp_solution(s) ....: for s in ....: chain( ....: solve_gevp_zero(GG, I, J), ....: solve_gevp_nonzero(GG, GG, GG, I, J) ....: ) ....: ) ....: == ....: set( ....: _normalize_gevp_solution(s) ....: for s in ....: _solve_gevp_naive(GG,GG,GG,I,J) ....: ) ....: for I in G_index_sets ....: for J in G_index_sets ....: ) True
>>> from sage.all import * >>> from itertools import chain >>> from sage.geometry.cone_critical_angles import ( ... _normalize_gevp_solution, ... _solve_gevp_naive, ... gevp_licis, ... solve_gevp_nonzero, ... solve_gevp_zero) >>> K = cones.schur(Integer(3)) >>> gs = [g.change_ring(AA).normalized() for g in K] >>> G = matrix.column(gs) >>> GG = G.transpose() * G >>> G_index_sets = list(gevp_licis(G)) >>> all( ... set( ... _normalize_gevp_solution(s) ... for s in ... chain( ... solve_gevp_zero(GG, I, J), ... solve_gevp_nonzero(GG, GG, GG, I, J) ... ) ... ) ... == ... set( ... _normalize_gevp_solution(s) ... for s in ... _solve_gevp_naive(GG,GG,GG,I,J) ... ) ... for I in G_index_sets ... for J in G_index_sets ... ) True
- sage.geometry.cone_critical_angles.solve_gevp_zero(M, I, J)[source]#
Solve the generalized eigenvalue problem in Theorem 3 [Or2020] for a zero eigenvalue using Propositions 3 and 4 [Or2020].
INPUT:
M
– the matrix whose \((i,j)\)-th entry is the inner product of \(g_{i}\) and \(h_{j}\) as in Proposition 6 [Or2020]I
– a linearly independent column-index set for the matrix \(G\) that appears in Theorem 3 [Or2020]J
– a linearly independent column-index set for the matrix \(H\) that appears in Theorem 3 [Or2020]
OUTPUT:
A generator of
(eigenvalue, xi, eta, multiplicity)
quartets whereeigenvalue
is zero (the eigenvalue of the system)xi
is the first (lengthlen(I)
) component of an eigenvector associated witheigenvalue
eta
is the second (lengthlen(J)
) component of an eigenvector associated witheigenvalue
multiplicity
is the dimension of the eigenspace associated witheigenvalue
ALGORITHM:
Proposition 4 in [Or2020] is used.
EXAMPLES:
This particular configuration results in the zero matrix in the eigenvalue problem, so the only solutions correspond to the eigenvalue zero:
sage: from sage.geometry.cone_critical_angles import solve_gevp_zero sage: K = cones.nonnegative_orthant(2) sage: G = matrix.column(K.rays()) sage: GG = G.transpose() * G sage: I = [0] sage: J = [1] sage: list(solve_gevp_zero(GG, I, J)) [(0, (1), (0), 2), (0, (0), (1), 2)]
>>> from sage.all import * >>> from sage.geometry.cone_critical_angles import solve_gevp_zero >>> K = cones.nonnegative_orthant(Integer(2)) >>> G = matrix.column(K.rays()) >>> GG = G.transpose() * G >>> I = [Integer(0)] >>> J = [Integer(1)] >>> list(solve_gevp_zero(GG, I, J)) [(0, (1), (0), 2), (0, (0), (1), 2)]