]>
gitweb.michael.orlitzky.com - dunshire.git/blob - src/dunshire/matrices.py
2 Utility functions for working with CVXOPT matrices (instances of the
3 class:`cvxopt.base.matrix` class).
7 from cvxopt
import matrix
8 from cvxopt
.lapack
import gees
, syevr
13 def append_col(left
, right
):
15 Append two matrices side-by-side.
21 The two matrices to append to one another.
27 A new matrix consisting of ``right`` appended to the right
33 >>> A = matrix([1,2,3,4], (2,2))
34 >>> B = matrix([5,6,7,8,9,10], (2,3))
43 >>> print(append_col(A,B))
49 return matrix([left
.trans(), right
.trans()]).trans()
52 def append_row(top
, bottom
):
54 Append two matrices top-to-bottom.
60 The two matrices to append to one another.
66 A new matrix consisting of ``bottom`` appended below ``top``.
71 >>> A = matrix([1,2,3,4], (2,2))
72 >>> B = matrix([5,6,7,8,9,10], (3,2))
82 >>> print(append_row(A,B))
91 return matrix([top
, bottom
])
94 def eigenvalues(symmat
):
96 Return the eigenvalues of the given symmetric real matrix.
98 On the surface, this appears redundant to the :func:`eigenvalues_re`
99 function. However, if we know in advance that our input is
100 symmetric, a better algorithm can be used.
106 The real symmetric matrix whose eigenvalues you want.
112 A list of the eigenvalues (in no particular order) of ``symmat``.
118 If the input matrix is not symmetric.
123 >>> A = matrix([[2,1],[1,2]], tc='d')
127 If the input matrix is not symmetric, it may not have real
128 eigenvalues, and we don't know what to do::
130 >>> A = matrix([[1,2],[3,4]])
132 Traceback (most recent call last):
134 TypeError: input must be a symmetric real matrix
137 if not norm(symmat
.trans() - symmat
) < options
.ABS_TOL
:
138 # Ensure that ``symmat`` is symmetric (and thus square).
139 raise TypeError('input must be a symmetric real matrix')
141 domain_dim
= symmat
.size
[0]
142 eigs
= matrix(0, (domain_dim
, 1), tc
='d')
147 def eigenvalues_re(anymat
):
149 Return the real parts of the eigenvalues of the given square matrix.
155 The square matrix whose eigenvalues you want.
161 A list of the real parts (in no particular order) of the
162 eigenvalues of ``anymat``.
168 If the input matrix is not square.
173 This is symmetric and has two real eigenvalues:
175 >>> A = matrix([[2,1],[1,2]], tc='d')
176 >>> sorted(eigenvalues_re(A))
179 But this rotation matrix has eigenvalues `i` and `-i`, both of whose
182 >>> A = matrix([[0,-1],[1,0]])
183 >>> eigenvalues_re(A)
186 If the input matrix is not square, it doesn't have eigenvalues::
188 >>> A = matrix([[1,2],[3,4],[5,6]])
189 >>> eigenvalues_re(A)
190 Traceback (most recent call last):
192 TypeError: input matrix must be square
195 if not anymat
.size
[0] == anymat
.size
[1]:
196 raise TypeError('input matrix must be square')
198 domain_dim
= anymat
.size
[0]
199 eigs
= matrix(0, (domain_dim
, 1), tc
='z')
201 # Create a copy of ``anymat`` here for two reasons:
203 # 1. ``gees`` clobbers its input.
204 # 2. We need to ensure that the type code of ``dummy`` is 'd' or 'z'.
206 dummy
= matrix(anymat
, anymat
.size
, tc
='d')
209 return [eig
.real
for eig
in eigs
]
212 def identity(domain_dim
):
214 Create an identity matrix of the given dimensions.
220 The dimension of the vector space on which the identity will act.
226 A ``domain_dim``-by-``domain_dim`` dense integer identity matrix.
232 If you ask for the identity on zero or fewer dimensions.
237 >>> print(identity(3))
245 raise ValueError('domain dimension must be positive')
247 entries
= [int(i
== j
)
248 for i
in range(domain_dim
)
249 for j
in range(domain_dim
)]
250 return matrix(entries
, (domain_dim
, domain_dim
))
253 def inner_product(vec1
, vec2
):
255 Compute the Euclidean inner product of two vectors.
261 The two vectors whose inner product you want.
267 The inner product of ``vec1`` and ``vec2``.
273 If the lengths of ``vec1`` and ``vec2`` differ.
280 >>> inner_product(x,y)
283 >>> x = matrix([1,1,1])
284 >>> y = matrix([2,3,4], (1,3))
285 >>> inner_product(x,y)
290 >>> inner_product(x,y)
291 Traceback (most recent call last):
293 TypeError: the lengths of vec1 and vec2 must match
296 if not len(vec1
) == len(vec2
):
297 raise TypeError('the lengths of vec1 and vec2 must match')
299 return sum([x
*y
for (x
, y
) in zip(vec1
, vec2
)])
302 def norm(matrix_or_vector
):
304 Return the Frobenius norm of a matrix or vector.
306 When the input is a vector, its matrix-Frobenius norm is the same
307 thing as its vector-Euclidean norm.
312 matrix_or_vector : matrix
313 The matrix or vector whose norm you want.
319 The norm of ``matrix_or_vector``.
324 >>> v = matrix([1,1])
325 >>> print('{:.5f}'.format(norm(v)))
328 >>> A = matrix([1,1,1,1], (2,2))
333 return sqrt(inner_product(matrix_or_vector
, matrix_or_vector
))
338 Create a long vector in column-major order from ``mat``.
344 Any sort of real matrix that you want written as a long vector.
350 An ``len(mat)``-by-``1`` long column vector containign the
351 entries of ``mat`` in column major order.
356 >>> A = matrix([[1,2],[3,4]])
369 Note that if ``mat`` is a vector, this function is a no-op:
371 >>> v = matrix([1,2,3,4], (4,1))
386 return matrix(mat
, (len(mat
), 1))