From b4b48ce895557cc5708b090568cfa33e7d09817f Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Fri, 30 Nov 2018 11:19:56 -0500 Subject: [PATCH] mjo/cone/symmetric_psd: don't use recusion; work on trivial matrices. --- mjo/cone/symmetric_psd.py | 49 +++++++++++++++++++++++++-------------- 1 file changed, 32 insertions(+), 17 deletions(-) diff --git a/mjo/cone/symmetric_psd.py b/mjo/cone/symmetric_psd.py index e75c152..40fd90c 100644 --- a/mjo/cone/symmetric_psd.py +++ b/mjo/cone/symmetric_psd.py @@ -46,6 +46,13 @@ def is_symmetric_psd(A): sage: is_symmetric_psd(A) False + The trivial matrix in a trivial space is trivially symmetric and + positive-semidefinite:: + + sage: A = matrix(QQ, 0,0) + sage: is_symmetric_psd(A) + True + """ if A.base_ring() == SR: @@ -60,7 +67,11 @@ def is_symmetric_psd(A): # semidefinite. For that we can consult its minimum eigenvalue, # which should be zero or greater. Since ``A`` is symmetric, its # eigenvalues are guaranteed to be real. - return min(A.eigenvalues()) >= 0 + if A.is_zero(): + # A is trivial... so trivially positive-semudefinite. + return True + else: + return min(A.eigenvalues()) >= 0 def unit_eigenvectors(A): @@ -304,19 +315,23 @@ def random_symmetric_psd(V, accept_zero=True, rank=None): # Use the one the user gave us. rank_A = rank - # Begin with the zero matrix, and add projectors to it if we have any. - A = V.zero().column()*V.zero().row() - - # Careful, begin at idx=1 so that we only generate a projector - # when rank_A is greater than zero. - while A.rank() < rank_A: - v = V.random_element() - A += v.column()*v.row() - - if accept_zero or not A.is_zero(): - # We either don't care what ``A`` is, or it's non-zero, so - # just return it. - return A - else: - # Uh oh, we need to generate a new one. - return random_symmetric_psd(V, accept_zero, rank) + if n == 0 and not accept_zero: + # We're gonna loop forever trying to satisfy this... + raise ValueError('You must have accept_zero=True when V is trivial') + + # Loop until we find a suitable "A" that will then be returned. + while True: + # Begin with the zero matrix, and add projectors to it if we + # have any. + A = matrix.zero(V.base_ring(), n, n) + + # Careful, begin at idx=1 so that we only generate a projector + # when rank_A is greater than zero. + while A.rank() < rank_A: + v = V.random_element() + A += v.column()*v.row() + + if accept_zero or not A.is_zero(): + # We either don't care what ``A`` is, or it's non-zero, so + # just return it. + return A -- 2.44.2