]>
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).
8 from cvxopt
import matrix
9 from cvxopt
.lapack
import gees
, 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')
148 def eigenvalues_re(anymat
):
150 Return the real parts of the eigenvalues of the given square matrix.
156 The square matrix whose eigenvalues you want.
162 A list of the real parts (in no particular order) of the
163 eigenvalues of ``anymat``.
169 If the input matrix is not square.
174 This is symmetric and has two real eigenvalues:
176 >>> A = matrix([[2,1],[1,2]], tc='d')
177 >>> sorted(eigenvalues_re(A))
180 But this rotation matrix has eigenvalues `i` and `-i`, both of whose
183 >>> A = matrix([[0,-1],[1,0]])
184 >>> eigenvalues_re(A)
187 If the input matrix is not square, it doesn't have eigenvalues::
189 >>> A = matrix([[1,2],[3,4],[5,6]])
190 >>> eigenvalues_re(A)
191 Traceback (most recent call last):
193 TypeError: input matrix must be square
196 if not anymat
.size
[0] == anymat
.size
[1]:
197 raise TypeError('input matrix must be square')
199 domain_dim
= anymat
.size
[0]
200 eigs
= matrix(0, (domain_dim
, 1), tc
='z')
202 # Create a copy of ``anymat`` here for two reasons:
204 # 1. ``gees`` clobbers its input.
205 # 2. We need to ensure that the type code of ``dummy`` is 'd' or 'z'.
207 dummy
= matrix(anymat
, anymat
.size
, tc
='d')
210 return [eig
.real
for eig
in eigs
]
213 def identity(domain_dim
):
215 Create an identity matrix of the given dimensions.
221 The dimension of the vector space on which the identity will act.
227 A ``domain_dim``-by-``domain_dim`` dense integer identity matrix.
233 If you ask for the identity on zero or fewer dimensions.
238 >>> print(identity(3))
246 raise ValueError('domain dimension must be positive')
248 entries
= [int(i
== j
)
249 for i
in range(domain_dim
)
250 for j
in range(domain_dim
)]
251 return matrix(entries
, (domain_dim
, domain_dim
))
254 def inner_product(vec1
, vec2
):
256 Compute the Euclidean inner product of two vectors.
262 The two vectors whose inner product you want.
268 The inner product of ``vec1`` and ``vec2``.
274 If the lengths of ``vec1`` and ``vec2`` differ.
281 >>> inner_product(x,y)
284 >>> x = matrix([1,1,1])
285 >>> y = matrix([2,3,4], (1,3))
286 >>> inner_product(x,y)
291 >>> inner_product(x,y)
292 Traceback (most recent call last):
294 TypeError: the lengths of vec1 and vec2 must match
297 if not len(vec1
) == len(vec2
):
298 raise TypeError('the lengths of vec1 and vec2 must match')
300 return sum([x
*y
for (x
, y
) in zip(vec1
, vec2
)])
303 def norm(matrix_or_vector
):
305 Return the Frobenius norm of a matrix or vector.
307 When the input is a vector, its matrix-Frobenius norm is the same
308 thing as its vector-Euclidean norm.
313 matrix_or_vector : matrix
314 The matrix or vector whose norm you want.
320 The norm of ``matrix_or_vector``.
325 >>> v = matrix([1,1])
326 >>> print('{:.5f}'.format(norm(v)))
329 >>> A = matrix([1,1,1,1], (2,2))
334 return sqrt(inner_product(matrix_or_vector
, matrix_or_vector
))
339 Create a long vector in column-major order from ``mat``.
345 Any sort of real matrix that you want written as a long vector.
351 An ``len(mat)``-by-``1`` long column vector containign the
352 entries of ``mat`` in column major order.
357 >>> A = matrix([[1,2],[3,4]])
370 Note that if ``mat`` is a vector, this function is a no-op:
372 >>> v = matrix([1,2,3,4], (4,1))
387 return matrix(mat
, (len(mat
), 1))