]>
gitweb.michael.orlitzky.com - dunshire.git/blob - matrices.py
6d5bbc8f667dce05feb93ec9fabcf42340f1157e
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 two matrices to append to one another.
28 A new matrix consisting of ``right`` appended to the right
34 >>> A = matrix([1,2,3,4], (2,2))
35 >>> B = matrix([5,6,7,8,9,10], (2,3))
44 >>> print(append_col(A,B))
50 return matrix([left
.trans(), right
.trans()]).trans()
53 def append_row(top
, bottom
):
55 Append two matrices top-to-bottom.
61 The two matrices to append to one another.
67 A new matrix consisting of ``bottom`` appended below ``top``.
72 >>> A = matrix([1,2,3,4], (2,2))
73 >>> B = matrix([5,6,7,8,9,10], (3,2))
83 >>> print(append_row(A,B))
92 return matrix([top
, bottom
])
95 def eigenvalues(symmat
):
97 Return the eigenvalues of the given symmetric real matrix.
99 On the surface, this appears redundant to the :func:`eigenvalues_re`
100 function. However, if we know in advance that our input is
101 symmetric, a better algorithm can be used.
107 The real symmetric matrix whose eigenvalues you want.
113 A list of the eigenvalues (in no particular order) of ``symmat``.
119 If the input matrix is not symmetric.
124 >>> A = matrix([[2,1],[1,2]], tc='d')
128 If the input matrix is not symmetric, it may not have real
129 eigenvalues, and we don't know what to do::
131 >>> A = matrix([[1,2],[3,4]])
133 Traceback (most recent call last):
135 TypeError: input must be a symmetric real matrix
138 if not norm(symmat
.trans() - symmat
) < options
.ABS_TOL
:
139 # Ensure that ``symmat`` is symmetric (and thus square).
140 raise TypeError('input must be a symmetric real matrix')
142 domain_dim
= symmat
.size
[0]
143 eigs
= matrix(0, (domain_dim
, 1), tc
='d')
145 # Create a copy of ``symmat`` here because ``syevr`` clobbers it.
151 def eigenvalues_re(anymat
):
153 Return the real parts of the eigenvalues of the given square matrix.
159 The square matrix whose eigenvalues you want.
165 A list of the real parts (in no particular order) of the
166 eigenvalues of ``anymat``.
172 If the input matrix is not square.
177 This is symmetric and has two real eigenvalues:
179 >>> A = matrix([[2,1],[1,2]], tc='d')
180 >>> sorted(eigenvalues_re(A))
183 But this rotation matrix has eigenvalues `i` and `-i`, both of whose
186 >>> A = matrix([[0,-1],[1,0]])
187 >>> eigenvalues_re(A)
190 If the input matrix is not square, it doesn't have eigenvalues::
192 >>> A = matrix([[1,2],[3,4],[5,6]])
193 >>> eigenvalues_re(A)
194 Traceback (most recent call last):
196 TypeError: input matrix must be square
199 if not anymat
.size
[0] == anymat
.size
[1]:
200 raise TypeError('input matrix must be square')
202 domain_dim
= anymat
.size
[0]
203 eigs
= matrix(0, (domain_dim
, 1), tc
='z')
205 # Create a copy of ``anymat`` here for two reasons:
207 # 1. ``gees`` clobbers its input.
208 # 2. We need to ensure that the type code of ``dummy`` is 'd' or 'z'.
210 dummy
= matrix(anymat
, anymat
.size
, tc
='d')
213 return [eig
.real
for eig
in eigs
]
216 def identity(domain_dim
, typecode
='i'):
218 Create an identity matrix of the given dimensions.
224 The dimension of the vector space on which the identity will act.
226 typecode : {'i', 'd', 'z'}, optional
227 The type code for the returned matrix, defaults to 'i' for integers.
228 Can also be 'd' for real double, or 'z' for complex double.
234 A ``domain_dim``-by-``domain_dim`` dense integer identity matrix.
240 If you ask for the identity on zero or fewer dimensions.
245 >>> print(identity(3))
253 raise ValueError('domain dimension must be positive')
255 entries
= [int(i
== j
)
256 for i
in range(domain_dim
)
257 for j
in range(domain_dim
)]
258 return matrix(entries
, (domain_dim
, domain_dim
), tc
=typecode
)
261 def inner_product(vec1
, vec2
):
263 Compute the Euclidean inner product of two vectors.
269 The two vectors whose inner product you want.
275 The inner product of ``vec1`` and ``vec2``.
281 If the lengths of ``vec1`` and ``vec2`` differ.
288 >>> inner_product(x,y)
291 >>> x = matrix([1,1,1])
292 >>> y = matrix([2,3,4], (1,3))
293 >>> inner_product(x,y)
298 >>> inner_product(x,y)
299 Traceback (most recent call last):
301 TypeError: the lengths of vec1 and vec2 must match
304 if not len(vec1
) == len(vec2
):
305 raise TypeError('the lengths of vec1 and vec2 must match')
307 return sum([x
*y
for (x
, y
) in zip(vec1
, vec2
)])
310 def norm(matrix_or_vector
):
312 Return the Frobenius norm of a matrix or vector.
314 When the input is a vector, its matrix-Frobenius norm is the same
315 thing as its vector-Euclidean norm.
320 matrix_or_vector : matrix
321 The matrix or vector whose norm you want.
327 The norm of ``matrix_or_vector``.
332 >>> v = matrix([1,1])
336 >>> A = matrix([1,1,1,1], (2,2))
341 return sqrt(inner_product(matrix_or_vector
, matrix_or_vector
))
346 Create a long vector in column-major order from ``mat``.
352 Any sort of real matrix that you want written as a long vector.
358 An ``len(mat)``-by-``1`` long column vector containign the
359 entries of ``mat`` in column major order.
364 >>> A = matrix([[1,2],[3,4]])
377 Note that if ``mat`` is a vector, this function is a no-op:
379 >>> v = matrix([1,2,3,4], (4,1))
394 return matrix(mat
, (len(mat
), 1))
397 def condition_number(mat
):
399 Return the condition number of the given matrix.
401 The condition number of a matrix quantifies how hard it is to do
402 numerical computation with that matrix. It is usually defined as
403 the ratio of the norm of the matrix to the norm of its inverse, and
404 therefore depends on the norm used. One way to compute the condition
405 number with respect to the 2-norm is as the ratio of the matrix's
406 largest and smallest singular values. Since we have easy access to
407 those singular values, that is the algorithm we use.
409 The larger the condition number is, the worse the matrix is.
414 The matrix whose condition number you want.
420 The nonnegative condition number of ``mat``.
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]