]>
gitweb.michael.orlitzky.com - dunshire.git/blob - dunshire/matrices.py
f35827d16be600cec3ee5d1c2d0fbcea27a9484d
2 Utility functions for working with CVXOPT matrices (instances of the
3 class:`cvxopt.base.matrix` class).
8 from cvxopt
import matrix
9 from cvxopt
.lapack
import gees
, gesdd
, syevr
14 def append_col(left
, right
):
16 Append two matrices side-by-side.
22 The "original" matrix, the one that will wind up on the left.
25 The matrix to be appended on the right of ``left``.
31 A new matrix consisting of ``right`` appended to the right
37 >>> A = matrix([1,2,3,4], (2,2))
38 >>> B = matrix([5,6,7,8,9,10], (2,3))
47 >>> print(append_col(A,B))
53 return matrix([left
.trans(), right
.trans()]).trans()
56 def append_row(top
, bottom
):
58 Append two matrices top-to-bottom.
64 The "original" matrix, the one that will wind up on top.
67 The matrix to be appended below ``top``.
73 A new matrix consisting of ``bottom`` appended below ``top``.
78 >>> A = matrix([1,2,3,4], (2,2))
79 >>> B = matrix([5,6,7,8,9,10], (3,2))
89 >>> print(append_row(A,B))
98 return matrix([top
, bottom
])
101 def eigenvalues(symmat
):
103 Return the eigenvalues of the given symmetric real matrix.
105 On the surface, this appears redundant to the :func:`eigenvalues_re`
106 function. However, if we know in advance that our input is
107 symmetric, a better algorithm can be used.
113 The real symmetric matrix whose eigenvalues you want.
119 A list of the eigenvalues (in no particular order) of ``symmat``.
125 If the input matrix is not symmetric.
130 >>> A = matrix([[2,1],[1,2]], tc='d')
134 If the input matrix is not symmetric, it may not have real
135 eigenvalues, and we don't know what to do::
137 >>> A = matrix([[1,2],[3,4]])
139 Traceback (most recent call last):
141 TypeError: input must be a symmetric real matrix
144 if not norm(symmat
.trans() - symmat
) < options
.ABS_TOL
:
145 # Ensure that ``symmat`` is symmetric (and thus square).
146 raise TypeError('input must be a symmetric real matrix')
148 domain_dim
= symmat
.size
[0]
149 eigs
= matrix(0, (domain_dim
, 1), tc
='d')
151 # Create a copy of ``symmat`` here because ``syevr`` clobbers it.
157 def eigenvalues_re(anymat
):
159 Return the real parts of the eigenvalues of the given square matrix.
165 The square matrix whose eigenvalues you want.
171 A list of the real parts (in no particular order) of the
172 eigenvalues of ``anymat``.
178 If the input matrix is not square.
183 This is symmetric and has two real eigenvalues:
185 >>> A = matrix([[2,1],[1,2]], tc='d')
186 >>> sorted(eigenvalues_re(A))
189 But this rotation matrix has eigenvalues `i` and `-i`, both of whose
192 >>> A = matrix([[0,-1],[1,0]])
193 >>> eigenvalues_re(A)
196 If the input matrix is not square, it doesn't have eigenvalues::
198 >>> A = matrix([[1,2],[3,4],[5,6]])
199 >>> eigenvalues_re(A)
200 Traceback (most recent call last):
202 TypeError: input matrix must be square
205 if not anymat
.size
[0] == anymat
.size
[1]:
206 raise TypeError('input matrix must be square')
208 domain_dim
= anymat
.size
[0]
209 eigs
= matrix(0, (domain_dim
, 1), tc
='z')
211 # Create a copy of ``anymat`` here for two reasons:
213 # 1. ``gees`` clobbers its input.
214 # 2. We need to ensure that the type code of ``dummy`` is 'd' or 'z'.
216 dummy
= matrix(anymat
, anymat
.size
, tc
='d')
219 return [eig
.real
for eig
in eigs
]
222 def identity(domain_dim
, typecode
='i'):
224 Create an identity matrix of the given dimensions.
230 The dimension of the vector space on which the identity will act.
232 typecode : {'i', 'd', 'z'}, optional
233 The type code for the returned matrix, defaults to 'i' for integers.
234 Can also be 'd' for real double, or 'z' for complex double.
240 A ``domain_dim``-by-``domain_dim`` dense integer identity matrix.
246 If you ask for the identity on zero or fewer dimensions.
251 >>> print(identity(3))
259 raise ValueError('domain dimension must be positive')
261 entries
= [int(i
== j
)
262 for i
in range(domain_dim
)
263 for j
in range(domain_dim
)]
264 return matrix(entries
, (domain_dim
, domain_dim
), tc
=typecode
)
267 def inner_product(vec1
, vec2
):
269 Compute the Euclidean inner product of two vectors.
275 The first vector, whose inner product with ``vec2`` you want.
278 The second vector, whose inner product with ``vec1`` you want.
284 The inner product of ``vec1`` and ``vec2``.
290 If the lengths of ``vec1`` and ``vec2`` differ.
297 >>> inner_product(x,y)
300 >>> x = matrix([1,1,1])
301 >>> y = matrix([2,3,4], (1,3))
302 >>> inner_product(x,y)
307 >>> inner_product(x,y)
308 Traceback (most recent call last):
310 TypeError: the lengths of vec1 and vec2 must match
313 if not len(vec1
) == len(vec2
):
314 raise TypeError('the lengths of vec1 and vec2 must match')
316 return sum([x
*y
for (x
, y
) in zip(vec1
, vec2
)])
319 def norm(matrix_or_vector
):
321 Return the Frobenius norm of a matrix or vector.
323 When the input is a vector, its matrix-Frobenius norm is the same
324 thing as its vector-Euclidean norm.
329 matrix_or_vector : matrix
330 The matrix or vector whose norm you want.
336 The norm of ``matrix_or_vector``.
341 >>> v = matrix([1,1])
345 >>> A = matrix([1,1,1,1], (2,2))
350 return sqrt(inner_product(matrix_or_vector
, matrix_or_vector
))
355 Return the spectral norm of a matrix.
357 The spectral norm of a matrix is its largest singular value, and it
358 corresponds to the operator norm induced by the vector Euclidean norm.
364 The matrix whose spectral norm you want.
368 >>> specnorm(identity(3))
371 >>> specnorm(5*identity(4))
375 num_eigs
= min(mat
.size
)
376 eigs
= matrix(0, (num_eigs
, 1), tc
='d')
378 if any([isinstance(entry
, complex) for entry
in mat
]):
380 dummy
= matrix(mat
, mat
.size
, tc
=typecode
)
391 Create a long vector in column-major order from ``mat``.
397 Any sort of real matrix that you want written as a long vector.
403 An ``len(mat)``-by-``1`` long column vector containign the
404 entries of ``mat`` in column major order.
409 >>> A = matrix([[1,2],[3,4]])
422 Note that if ``mat`` is a vector, this function is a no-op:
424 >>> v = matrix([1,2,3,4], (4,1))
439 return matrix(mat
, (len(mat
), 1))
442 def condition_number(mat
):
444 Return the condition number of the given matrix.
446 The condition number of a matrix quantifies how hard it is to do
447 numerical computation with that matrix. It is usually defined as
448 the ratio of the norm of the matrix to the norm of its inverse, and
449 therefore depends on the norm used. One way to compute the condition
450 number with respect to the 2-norm is as the ratio of the matrix's
451 largest and smallest singular values. Since we have easy access to
452 those singular values, that is the algorithm we use.
454 The larger the condition number is, the worse the matrix is.
459 The matrix whose condition number you want.
465 The nonnegative condition number of ``mat``.
470 >>> condition_number(identity(3))
473 >>> A = matrix([[2,1],[1,2]])
474 >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
477 >>> A = matrix([[2,1j],[-1j,2]])
478 >>> abs(condition_number(A) - 3.0) < options.ABS_TOL
482 num_eigs
= min(mat
.size
)
483 eigs
= matrix(0, (num_eigs
, 1), tc
='d')
485 if any([isinstance(entry
, complex) for entry
in mat
]):
487 dummy
= matrix(mat
, mat
.size
, tc
=typecode
)
491 return eigs
[0]/eigs
[-1]