]>
gitweb.michael.orlitzky.com - dunshire.git/blob - 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
, gesdd
, 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
, typecode
='i'):
214 Create an identity matrix of the given dimensions.
220 The dimension of the vector space on which the identity will act.
222 typecode : {'i', 'd', 'z'}, optional
223 The type code for the returned matrix, defaults to 'i' for integers.
224 Can also be 'd' for real double, or 'z' for complex double.
230 A ``domain_dim``-by-``domain_dim`` dense integer identity matrix.
236 If you ask for the identity on zero or fewer dimensions.
241 >>> print(identity(3))
249 raise ValueError('domain dimension must be positive')
251 entries
= [int(i
== j
)
252 for i
in range(domain_dim
)
253 for j
in range(domain_dim
)]
254 return matrix(entries
, (domain_dim
, domain_dim
), tc
=typecode
)
257 def inner_product(vec1
, vec2
):
259 Compute the Euclidean inner product of two vectors.
265 The two vectors whose inner product you want.
271 The inner product of ``vec1`` and ``vec2``.
277 If the lengths of ``vec1`` and ``vec2`` differ.
284 >>> inner_product(x,y)
287 >>> x = matrix([1,1,1])
288 >>> y = matrix([2,3,4], (1,3))
289 >>> inner_product(x,y)
294 >>> inner_product(x,y)
295 Traceback (most recent call last):
297 TypeError: the lengths of vec1 and vec2 must match
300 if not len(vec1
) == len(vec2
):
301 raise TypeError('the lengths of vec1 and vec2 must match')
303 return sum([x
*y
for (x
, y
) in zip(vec1
, vec2
)])
306 def norm(matrix_or_vector
):
308 Return the Frobenius norm of a matrix or vector.
310 When the input is a vector, its matrix-Frobenius norm is the same
311 thing as its vector-Euclidean norm.
316 matrix_or_vector : matrix
317 The matrix or vector whose norm you want.
323 The norm of ``matrix_or_vector``.
328 >>> v = matrix([1,1])
329 >>> print('{:.5f}'.format(norm(v)))
332 >>> A = matrix([1,1,1,1], (2,2))
337 return sqrt(inner_product(matrix_or_vector
, matrix_or_vector
))
342 Create a long vector in column-major order from ``mat``.
348 Any sort of real matrix that you want written as a long vector.
354 An ``len(mat)``-by-``1`` long column vector containign the
355 entries of ``mat`` in column major order.
360 >>> A = matrix([[1,2],[3,4]])
373 Note that if ``mat`` is a vector, this function is a no-op:
375 >>> v = matrix([1,2,3,4], (4,1))
390 return matrix(mat
, (len(mat
), 1))
393 def condition_number(mat
):
395 Return the condition number of the given matrix.
397 The condition number of a matrix quantifies how hard it is to do
398 numerical computation with that matrix. It is usually defined as
399 the ratio of the norm of the matrix to the norm of its inverse, and
400 therefore depends on the norm used. One way to compute the condition
401 number with respect to the 2-norm is as the ratio of the matrix's
402 largest and smallest singular values. Since we have easy access to
403 those singular values, that is the algorithm we use.
405 The larger the condition number is, the worse the matrix is.
410 The matrix whose condition number you want.
416 The nonnegative condition number of ``mat``.
421 >>> condition_number(identity(1, typecode='d'))
423 >>> condition_number(identity(2, typecode='d'))
425 >>> condition_number(identity(3, typecode='d'))
428 >>> A = matrix([[2,1],[1,2]], tc='d')
429 >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
432 >>> A = matrix([[2,1j],[-1j,2]], tc='z')
433 >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
437 num_eigs
= min(mat
.size
)
438 eigs
= matrix(0, (num_eigs
,1), tc
='d')
442 return eigs
[0]/eigs
[-1]