CVXOPT User’s Guide¶

Dense Matrices¶

A dense matrix is created by calling the function matrix. The arguments specify the values of the coefficients, the dimensions, and the type (integer, double, or complex) of the matrix.

cvxopt.matrix(x[, size[, tc]])¶

size is a tuple of length two with the matrix dimensions. The number of rows and/or the number of columns can be zero.

tc stands for type code. The possible values are 'i', 'd', and 'z', for integer, real (double), and complex matrices, respectively.

x can be a number, a sequence of numbers, a dense or sparse matrix, a one- or two-dimensional NumPy array, or a list of lists of matrices and numbers.

• If x is a number (Python integer, float, or complex number), a matrix is created with the dimensions specified by size and with all the coefficients equal to x. The default value of size is (1,1), and the default value of tc is the type of x. If necessary, the type of x is converted (from integer to double when used to create a matrix of type 'd', and from integer or double to complex when used to create a matrix of type 'z').

  >>> from cvxopt import matrix
  >>> A = matrix(1, (1,4))
  >>> print(A)
  [ 1  1  1  1]
  >>> A = matrix(1.0, (1,4))
  >>> print(A)
  [ 1.00e+00  1.00e+00  1.00e+00  1.00e+00]
  >>> A = matrix(1+1j)
  >>> print(A)
  [ 1.00e+00+j1.00e+00]
  

• If x is a sequence of numbers (list, tuple, array, array array, …), then the numbers are interpreted as the coefficients of a matrix in column-major order. The length of x must be equal to the product of size[0] and size[1]. If size is not specified, a matrix with one column is created. If tc is not specified, it is determined from the elements of x (and if that is impossible, for example because x is an empty list, a value 'i' is used). Type conversion takes place as for scalar x.

  The following example shows several ways to define the same integer matrix.

  >>> A = matrix([0, 1, 2, 3], (2,2))
  >>> A = matrix((0, 1, 2, 3), (2,2))
  >>> A = matrix(range(4), (2,2))
  >>> from array import array
  >>> A = matrix(array('i', [0,1,2,3]), (2,2))
  >>> print(A)
  [ 0  2]
  [ 1  3]
  

  In Python 2.7 the following also works.

  >>> A = matrix(xrange(4), (2,2))
  

• If x is a dense or sparse matrix, then the coefficients of x are copied, in column-major order, to a new matrix of the given size. The total number of elements in the new matrix (the product of size[0] and size[1]) must be the same as the product of the dimensions of x. If size is not specified, the dimensions of x are used. The default value of tc is the type of x. Type conversion takes place when the type of x differs from tc, in a similar way as for scalar x.

  >>> A = matrix([1., 2., 3., 4., 5., 6.], (2,3))
  >>> print(A)
  [ 1.00e+00  3.00e+00  5.00e+00]
  [ 2.00e+00  4.00e+00  6.00e+00]
  >>> B = matrix(A, (3,2))
  >>> print(B)
  [ 1.00e+00  4.00e+00]
  [ 2.00e+00  5.00e+00]
  [ 3.00e+00  6.00e+00]
  >>> C = matrix(B, tc='z')
  >>> print(C)
  [ 1.00e+00-j0.00e+00  4.00e+00-j0.00e+00]
  [ 2.00e+00-j0.00e+00  5.00e+00-j0.00e+00]
  [ 3.00e+00-j0.00e+00  6.00e+00-j0.00e+00]
  

  NumPy arrays can be converted to matrices.

  >>> from numpy import array
  >>> x = array([[1., 2., 3.], [4., 5., 6.]])
  >>> x
  array([[ 1.  2.  3.]
         [ 4.  5.  6.]])
  >>> print(matrix(x))
  [ 1.00e+00  2.00e+00  3.00e+00]
  [ 4.00e+00  5.00e+00  6.00e+00]
  

• If x is a list of lists of dense or sparse matrices and numbers (Python integer, float, or complex), then each element of x is interpreted as a block-column stored in column-major order. If size is not specified, the block-columns are juxtaposed to obtain a matrix with len(x) block-columns. If size is specified, then the matrix with len(x) block-columns is resized by copying its elements in column-major order into a matrix of the dimensions given by size. If tc is not specified, it is determined from the elements of x (and if that is impossible, for example because x is a list of empty lists, a value 'i' is used). The same rules for type conversion apply as for scalar x.

  >>> print(matrix([[1., 2.], [3., 4.], [5., 6.]]))
  [ 1.00e+00  3.00e+00  5.00e+00]
  [ 2.00e+00  4.00e+00  6.00e+00]
  >>> A1 = matrix([1, 2], (2,1))
  >>> B1 = matrix([6, 7, 8, 9, 10, 11], (2,3))
  >>> B2 = matrix([12, 13, 14, 15, 16, 17], (2,3))
  >>> B3 = matrix([18, 19, 20], (1,3))
  >>> C = matrix([[A1, 3.0, 4.0, 5.0], [B1, B2, B3]])
  >>> print(C)
  [ 1.00e+00  6.00e+00  8.00e+00  1.00e+01]
  [ 2.00e+00  7.00e+00  9.00e+00  1.10e+01]
  [ 3.00e+00  1.20e+01  1.40e+01  1.60e+01]
  [ 4.00e+00  1.30e+01  1.50e+01  1.70e+01]
  [ 5.00e+00  1.80e+01  1.90e+01  2.00e+01]
  

  A matrix with a single block-column can be represented by a single list (i.e., if x is a list of lists, and has length one, then the argument x can be replaced by x[0]).

  >>> D = matrix([B1, B2, B3])
  >>> print(D)
  [  6   8  10]
  [  7   9  11]
  [ 12  14  16]
  [ 13  15  17]
  [ 18  19  20]
  

Sparse Matrices¶

A general spmatrix object can be thought of as a triplet description of a sparse matrix, i.e., a list of entries of the matrix, with for each entry the value, row index, and column index. Entries that are not included in the list are assumed to be zero. For example, the sparse matrix

has the triplet description

The list may include entries with a zero value, so triplet descriptions are not necessarily unique. The list

is another triplet description of the same matrix.

An spmatrix object corresponds to a particular triplet description of a sparse matrix. We will refer to the entries in the triplet description as the nonzero entries of the object, even though they may have a numerical value zero.

Three functions are provided to create sparse matrices. The first, spmatrix, constructs a sparse matrix from a triplet description.

cvxopt.spmatrix(x, I, J[, size[, tc]])¶

I and J are sequences of integers (lists, tuples, array arrays, …) or integer matrices (matrix objects with typecode 'i'), containing the row and column indices of the nonzero entries. The lengths of I and J must be equal. If they are matrices, they are treated as lists of indices stored in column-major order, i.e., as lists list(I), respectively, list(J).

size is a tuple of nonnegative integers with the row and column dimensions of the matrix. The size argument is only needed when creating a matrix with a zero last row or last column. If size is not specified, it is determined from I and J: the default value for size[0] is max(I)+1 if I is nonempty and zero otherwise. The default value for size[1] is max(J)+1 if J is nonempty and zero otherwise.

tc is the typecode, 'd' or 'z', for double and complex matrices, respectively. Integer sparse matrices are not implemented.

x can be a number, a sequence of numbers, or a dense matrix. This argument specifies the numerical values of the nonzero entries.

• If x is a number (Python integer, float, or complex), a matrix is created with the sparsity pattern defined by I and J, and nonzero entries initialized to the value of x. The default value of tc is 'd' if x is integer or float, and 'z' if x is complex.

  The following code creates a 4 by 4 sparse identity matrix.

  >>> from cvxopt import spmatrix
  >>> A = spmatrix(1.0, range(4), range(4))
  >>> print(A)
     [ 1.00e+00     0         0         0    ]
     [    0      1.00e+00     0         0    ]
     [    0         0      1.00e+00     0    ]
     [    0         0         0      1.00e+00]
  

• If x is a sequence of numbers, a sparse matrix is created with the entries of x copied to the entries indexed by I and J. The list x must have the same length as I and J. The default value of tc is determined from the elements of x: 'd' if x contains integers and floating-point numbers or if x is an empty list, and 'z' if x contains at least one complex number.

  >>> A = spmatrix([2,-1,2,-2,1,4,3], [1,2,0,2,3,2,0], [0,0,1,1,2,3,4])
  >>> print(A)
    [    0      2.00e+00     0         0      3.00e+00]
    [ 2.00e+00     0         0         0         0    ]
    [-1.00e+00 -2.00e+00     0      4.00e+00     0    ]
    [    0         0      1.00e+00     0         0    ]
  

• If x is a dense matrix, a sparse matrix is created with all the entries of x copied, in column-major order, to the entries indexed by I and J. The matrix x must have the same length as I and J. The default value of tc is 'd' if x is an 'i' or 'd' matrix, and 'z' otherwise. If I and J contain repeated entries, the corresponding values of the coefficients are added.

The function sparse constructs a sparse matrix from a block-matrix description.

cvxopt.sparse(x[, tc])¶

tc is the typecode, 'd' or 'z', for double and complex matrices, respectively.

x can be a matrix, spmatrix, or a list of lists of matrices (matrix or spmatrix objects) and numbers (Python integer, float, or complex).

• If x is a matrix or spmatrix object, then a sparse matrix of the same size and the same numerical value is created. Numerical zeros in x are treated as structural zeros and removed from the triplet description of the new sparse matrix.
• If x is a list of lists of matrices (matrix or spmatrix objects) and numbers (Python integer, float, or complex) then each element of x is interpreted as a (block-)column matrix stored in colum-major order, and a block-matrix is constructed by juxtaposing the len(x) block-columns (as in matrix). Numerical zeros are removed from the triplet description of the new matrix.

>>> from cvxopt import matrix, spmatrix, sparse
>>> A = matrix([[1., 2., 0.], [2., 1., 2.], [0., 2., 1.]])
>>> print(A)
[ 1.00e+00  2.00e+00  0.00e+00]
[ 2.00e+00  1.00e+00  2.00e+00]
[ 0.00e+00  2.00e+00  1.00e+00]
>>> B = spmatrix([], [], [], (3,3))
>>> print(B)
[0 0 0]
[0 0 0]
[0 0 0]
>>> C = spmatrix([3, 4, 5], [0, 1, 2], [0, 1, 2])
>>> print(C)
[ 3.00e+00     0         0    ]
[    0      4.00e+00     0    ]
[    0         0      5.00e+00]
>>> D = sparse([[A, B], [B, C]])
>>> print(D)
[ 1.00e+00  2.00e+00     0         0         0         0    ]
[ 2.00e+00  1.00e+00  2.00e+00     0         0         0    ]
[    0      2.00e+00  1.00e+00     0         0         0    ]
[    0         0         0      3.00e+00     0         0    ]
[    0         0         0         0      4.00e+00     0    ]
[    0         0         0         0         0      5.00e+00]

A matrix with a single block-column can be represented by a single list.

>>> D = sparse([A, C])
>>> print(D)
[ 1.00e+00  2.00e+00     0    ]
[ 2.00e+00  1.00e+00  2.00e+00]
[    0      2.00e+00  1.00e+00]
[ 3.00e+00     0         0    ]
[    0      4.00e+00     0    ]
[    0         0      5.00e+00]

The function spdiag constructs a block-diagonal sparse matrix from a list of matrices.

cvxopt.spdiag(x)¶

x is a dense or sparse matrix with a single row or column, or a list of square dense or sparse matrices or scalars. If x is a matrix, a sparse diagonal matrix is returned with the entries of x on its diagonal. If x is list, a sparse block-diagonal matrix is returned with the elements in the list as its diagonal blocks.

>>> from cvxopt import matrix, spmatrix, spdiag
>>> A = 3.0
>>> B = matrix([[1,-2],[-2,1]])
>>> C = spmatrix([1,1,1,1,1],[0,1,2,0,0,],[0,0,0,1,2])
>>> D = spdiag([A, B, C])
>>> print(D)
[ 3.00e+00     0         0         0         0         0    ]
[    0      1.00e+00 -2.00e+00     0         0         0    ]
[    0     -2.00e+00  1.00e+00     0         0         0    ]
[    0         0         0      1.00e+00  1.00e+00  1.00e+00]
[    0         0         0      1.00e+00     0         0    ]
[    0         0         0      1.00e+00     0         0    ]

Arithmetic Operations¶

The following table lists the arithmetic operations defined for dense and sparse matrices. In the table A and B are dense or sparse matrices of compatible dimensions, c is a scalar (a Python number or a dense 1 by 1 matrix), D is a dense matrix, and e is a Python number.

The type of the result of these operations generally follows the Python conventions. For example, if A and c are integer, then in Python 2 the division A/c is interpreted as integer division and results in a type 'i' matrix, while in Python 3 it is interpreted as standard divison and results in a type 'd' matrix. An exception to the Python conventions is elementwise exponentiation: if D is an integer matrix and e is an integer number than D**e is a matrix of type 'd'.

Addition, subtraction, and matrix multiplication with two matrix operands result in a sparse matrix if both matrices are sparse, and in a dense matrix otherwise. The result of a scalar multiplication or division is dense if A is dense, and sparse if A is sparse. Postmultiplying a matrix with a number c means the same as premultiplying, i.e., scalar multiplication. Dividing a matrix by c means dividing all its entries by c.

If c in the expressions A+c, c+A, A-c, c-A is a number, then it is interpreted as a dense matrix with the same dimensions as A, type given by the type of c, and all its entries equal to c. If c is a 1 by 1 dense matrix and A is not 1 by 1, then c is interpreted as a dense matrix with the same size of A and all entries equal to c[0].

If c is a 1 by 1 dense matrix, then, if possible, the products c*A and A*c are interpreted as matrix-matrix products. If the product cannot be interpreted as a matrix-matrix product because the dimensions of A are incompatible, then the product is interpreted as the scalar multiplication with c[0]. The division A/c and remainder A%c with c a 1 by 1 matrix are always interpreted as A/c[0], resp., A%c[0].

The following in-place operations are also defined, but only if they do not change the type (sparse or dense, integer, real, or complex) of the matrix A. These in-place operations do not return a new matrix but modify the existing object A.

For example, if A has typecode 'i', then A += B is allowed if B has typecode 'i'. It is not allowed if B has typecode 'd' or 'z' because the addition A+B results in a 'd' or 'z' matrix and therefore cannot be assigned to A without changing its type. As another example, if A is a sparse matrix, then A += 1.0 is not allowed because the operation A = A + 1.0 results in a dense matrix, so it cannot be assigned to A without changing its type.

In-place matrix-matrix products are not allowed. (Except when c is a 1 by 1 dense matrix, in which case A *= c is interpreted as the scalar product A *= c[0].)

In-place remainder is only defined for dense A.

It is important to know when a matrix operation creates a new object. The following rules apply.

• A simple assignment (A = B) is given the standard Python interpretation, i.e., it assigns to the variable A a reference (or pointer) to the object referenced by B.

  >>> B = matrix([[1.,2.], [3.,4.]])
  >>> print(B)
  [ 1.00e+00  3.00e+00]
  [ 2.00e+00  4.00e+00]
  >>> A = B
  >>> A[0,0] = -1
  >>> print(B)   # modifying A[0,0] also modified B[0,0]
  [-1.00e+00  3.00e+00]
  [ 2.00e+00  4.00e+00]
  

• The regular (i.e., not in-place) arithmetic operations always return new objects.

  >>> B = matrix([[1.,2.], [3.,4.]])
  >>> A = +B
  >>> A[0,0] = -1
  >>> print(B)   # modifying A[0,0] does not modify B[0,0]
  [ 1.00e+00  3.00e+00]
  [ 2.00e+00  4.00e+00]
  

• The in-place operations directly modify the coefficients of the existing matrix object and do not create a new object.

  >>> B = matrix([[1.,2.], [3.,4.]])
  >>> A = B
  >>> A *= 2
  >>> print(B)   # in-place operation also changed B
  [ 2.00e+00  6.00e+00]
  [ 4.00e+00  8.00e+00]
  >>> A = 2*A
  >>> print(B)   # regular operation creates a new A, so does not change B
  [ 2.00e+00  6.00e+00]
  [ 4.00e+00  8.00e+00]
  

Indexing and Slicing¶

Matrices can be indexed using one or two arguments. In single-argument indexing of a matrix A, the index runs from -len(A) to len(A)-1, and is interpreted as an index in the one-dimensional array of coefficients of A in column-major order. Negative indices have the standard Python interpretation: for negative k, A[k] is the same element as A[len(A) + k].

Four different types of one-argument indexing are implemented.

1. The index can be a single integer. This returns a number, e.g., A[0] is the first element of A.
2. The index can be an integer matrix. This returns a column matrix: the command A[matrix([0,1,2,3])] returns the 4 by 1 matrix consisting of the first four elements of A. The size of the index matrix is ignored: A[matrix([0,1,2,3], (2,2))] returns the same 4 by 1 matrix.
3. The index can be a list of integers. This returns a column matrix, e.g., A[[0,1,2,3]] is the 4 by 1 matrix consisting of elements 0, 1, 2, 3 of A.
4. The index can be a Python slice. This returns a matrix with one column (possibly 0 by 1, or 1 by 1). For example, A[::2] is the column matrix defined by taking every other element of A, stored in column-major order. A[0:0] is a matrix with size (0,1).

Thus, single-argument indexing returns a scalar (if the index is an integer), or a matrix with one column. This is consistent with the interpretation that single-argument indexing accesses the matrix in column-major order.

Note that an index list or an index matrix are equivalent, but they are both useful, especially when we perform operations on index sets. For example, if I and J are lists then I+J is the concatenated list, and 2*I is I repeated twice. If they are matrices, these operations are interpreted as arithmetic operations. For large index sets, indexing with integer matrices is also faster than indexing with lists.

The following example illustrates one-argument indexing.

>>> from cvxopt import matrix, spmatrix
>>> A = matrix(range(16), (4,4), 'd')
>>> print(A)
[ 0.00e+00  4.00e+00  8.00e+00  1.20e+01]
[ 1.00e+00  5.00e+00  9.00e+00  1.30e+01]
[ 2.00e+00  6.00e+00  1.00e+01  1.40e+01]
[ 3.00e+00  7.00e+00  1.10e+01  1.50e+01]
>>> A[4]
4.0
>>> I = matrix([0, 5, 10, 15])
>>> print(A[I])      # the diagonal
[ 0.00e+00]
[ 5.00e+00]
[ 1.00e+01]
[ 1.50e+01]
>>> I = [0,2];  J = [1,3]
>>> print(A[2*I+J])  # duplicate I and append J
[ 0.00e+00]
[ 2.00e+00]
[ 0.00e+00]
[ 2.00e+00]
[ 1.00e+00]
[ 3.00e+00]
>>> I = matrix([0, 2]);  J =  matrix([1, 3])
>>> print(A[2*I+J])  # multiply I by 2 and add J
[ 1.00e+00]
[ 7.00e+00]
>>> print(A[4::4])   # get every fourth element skipping the first four
[ 4.00e+00]
[ 8.00e+00]
[ 1.20e+01]

In two-argument indexing the arguments can be any combinations of the four types listed above. The first argument indexes the rows of the matrix and the second argument indexes the columns. If both indices are scalars, then a scalar is returned. In all other cases, a matrix is returned. We continue the example.

>>> print(A[:,1])
[ 4.00e+00]
[ 5.00e+00]
[ 6.00e+00]
[ 7.00e+00]
>>> J = matrix([0, 2])
>>> print(A[J,J])
[ 0.00e+00  8.00e+00]
[ 2.00e+00  1.00e+01]
>>> print(A[:2, -2:])
[ 8.00e+00  1.20e+01]
[ 9.00e+00  1.30e+01]
>>> A = spmatrix([0,2,-1,2,-2,1], [0,1,2,0,2,1], [0,0,0,1,1,2])
>>> print(A[:, [0,1]])
[ 0.00e+00  2.00e+00]
[ 2.00e+00     0    ]
[-1.00e+00 -2.00e+00]
>>> B = spmatrix([0,2*1j,0,-2], [1,2,1,2], [0,0,1,1,])
>>> print(B[-2:,-2:])
[ 0.00e+00-j0.00e+00  0.00e+00-j0.00e+00]
[ 0.00e+00+j2.00e+00 -2.00e+00-j0.00e+00]

Expressions of the form A[I] or A[I,J] can also appear on the left-hand side of an assignment. The right-hand side must be a scalar (i.e., a number or a 1 by 1 dense matrix), a sequence of numbers, or a dense or sparse matrix. If the right-hand side is a scalar, it is interpreted as a dense matrix with identical entries and the dimensions of the left-hand side. If the right-hand side is a sequence of numbers (list, tuple, array array, range object, …) its values are interpreted as the coefficients of a dense matrix in column-major order. If the right-hand side is a matrix (matrix or spmatrix), it must have the same size as the left-hand side. Sparse matrices are converted to dense in the assignment to a dense matrix.

Indexed assignments are only allowed if they do not change the type of the matrix. For example, if A is a matrix with type 'd', then A[I] = B is only permitted if B is an integer, a float, or a matrix of type 'i' or 'd'. If A is an integer matrix, then A[I] = B is only permitted if B is an integer or an integer matrix.

The following examples illustrate indexed assignment.

>>> A = matrix(range(16), (4,4))
>>> A[::2,::2] = matrix([[-1, -2], [-3, -4]])
>>> print(A)
[ -1   4  -3  12]
[  1   5   9  13]
[ -2   6  -4  14]
[  3   7  11  15]
>>> A[::5] += 1
>>> print(A)
[  0   4  -3  12]
[  1   6   9  13]
[ -2   6  -3  14]
[  3   7  11  16]
>>> A[0,:] = -1, 1, -1, 1
>>> print(A)
[ -1   1  -1   1]
[  1   6   9  13]
[ -2   6  -3  14]
[  3   7  11  16]
>>> A[2:,2:] = range(4)
>>> print(A)
[ -1   1  -1   1]
[  1   6   9  13]
[ -2   6   0   2]
[  3   7   1   3]
>>> A = spmatrix([0,2,-1,2,-2,1], [0,1,2,0,2,1], [0,0,0,1,1,2])
>>> print(A)
[ 0.00e+00  2.00e+00     0    ]
[ 2.00e+00     0      1.00e+00]
[-1.00e+00 -2.00e+00     0    ]
>>> C = spmatrix([10,-20,30], [0,2,1], [0,0,1])
>>> print(C)
[ 1.00e+01     0    ]
[    0      3.00e+01]
[-2.00e+01     0    ]
>>> A[:,0] = C[:,0]
>>> print(A)
[ 1.00e+01  2.00e+00     0    ]
[    0         0      1.00e+00]
[-2.00e+01 -2.00e+00     0    ]
>>> D = matrix(range(6), (3,2))
>>> A[:,0] = D[:,0]
>>> print(A)
[ 0.00e+00  2.00e+00     0    ]
[ 1.00e+00     0      1.00e+00]
[ 2.00e+00 -2.00e+00     0    ]
>>> A[:,0] = 1
>>> print(A)
[ 1.00e+00  2.00e+00     0    ]
[ 1.00e+00     0      1.00e+00]
[ 1.00e+00 -2.00e+00     0    ]
>>> A[:,0] = 0
>>> print(A)
[ 0.00e+00  2.00e+00     0    ]
[ 0.00e+00     0      1.00e+00]
[ 0.00e+00 -2.00e+00     0    ]

Attributes and Methods¶

Dense and sparse matrices have the following attributes.

size¶

A tuple with the dimensions of the matrix. The size of the matrix can be changed by altering this attribute, as long as the number of elements in the matrix remains unchanged.

typecode¶

A character, either 'i', 'd', or 'z', for integer, real, and complex matrices, respectively. A read-only attribute.

trans()¶

Returns the transpose of the matrix as a new matrix. One can also use A.T instead of A.trans().

ctrans()¶

Returns the conjugate transpose of the matrix as a new matrix. One can also use A.H instead of A.ctrans().

real()¶

For complex matrices, returns the real part as a real matrix. For integer and real matrices, returns a copy of the matrix.

imag()¶

For complex matrices, returns the imaginary part as a real matrix. For integer and real matrices, returns an integer or real zero matrix.

In addition, sparse matrices have the following attributes.

V¶

A single-column dense matrix containing the numerical values of the nonzero entries in column-major order. Making an assignment to the attribute is an efficient way of changing the values of the sparse matrix, without changing the sparsity pattern.

When the attribute V is read, a copy of V is returned, as a new dense matrix. This implies, for example, that an indexed assignment A.V[I] = B does not work, or at least cannot be used to modify A. Instead the attribute V will be read and returned as a new matrix; then the elements of this new matrix are modified.

I¶

A single-column integer dense matrix with the row indices of the entries in V. A read-only attribute.

J¶

A single-column integer dense matrix with the column indices of the entries in V. A read-only attribute.

CCS¶

A triplet (column pointers, row indices, values) with the compressed-column-storage representation of the matrix. A read-only attribute. This attribute can be used to export sparse matrices to other packages such as MOSEK.

The next example below illustrates assignments to V.

>>> from cvxopt import spmatrix, matrix
>>> A = spmatrix(range(5), [0,1,1,2,2], [0,0,1,1,2])
>>> print(A)
[ 0.00e+00     0         0    ]
[ 1.00e+00  2.00e+00     0    ]
[    0      3.00e+00  4.00e+00]
>>> B = spmatrix(A.V, A.J, A.I, (4,4))  # transpose and add a zero row and column
>>> print(B)
[ 0.00e+00  1.00e+00     0         0    ]
[    0      2.00e+00  3.00e+00     0    ]
[    0         0      4.00e+00     0    ]
[    0         0         0         0    ]
>>> B.V = matrix([1., 7., 8., 6., 4.])   # assign new values to nonzero entries
>>> print(B)
[ 1.00e+00  7.00e+00     0         0    ]
[    0      8.00e+00  6.00e+00     0    ]
[    0         0      4.00e+00     0    ]
[    0         0         0         0    ]

The following attributes and methods are defined for dense matrices.

tofile(f)¶

Writes the elements of the matrix in column-major order to a binary file f.

fromfile(f)¶

Reads the contents of a binary file f into the matrix object.

The last two methods are illustrated in the following examples.

>>> from cvxopt import matrix, spmatrix
>>> A = matrix([[1.,2.,3.], [4.,5.,6.]])
>>> print(A)
[ 1.00e+00  4.00e+00]
[ 2.00e+00  5.00e+00]
[ 3.00e+00  6.00e+00]
>>> f = open('mat.bin','wb')
>>> A.tofile(f)
>>> f.close()
>>> B = matrix(0.0, (2,3))
>>> f = open('mat.bin','rb')
>>> B.fromfile(f)
>>> f.close()
>>> print(B)
[ 1.00e+00  3.00e+00  5.00e+00]
[ 2.00e+00  4.00e+00  6.00e+00]
>>> A = spmatrix(range(5), [0,1,1,2,2], [0,0,1,1,2])
>>> f = open('test.bin','wb')
>>> A.V.tofile(f)
>>> A.I.tofile(f)
>>> A.J.tofile(f)
>>> f.close()
>>> f = open('test.bin','rb')
>>> V = matrix(0.0, (5,1));  V.fromfile(f)
>>> I = matrix(0, (5,1));  I.fromfile(f)
>>> J = matrix(0, (5,1));  J.fromfile(f)
>>> B = spmatrix(V, I, J)
>>> print(B)
[ 0.00e+00     0         0    ]
[ 1.00e+00  2.00e+00     0    ]
[    0      3.00e+00  4.00e+00]

Note that the dump and load functions in the pickle module offer a convenient alternative for writing matrices to files and reading matrices from files.

Built-In Functions¶

Many Python built-in functions and operations can be used with matrix arguments. We list some useful examples.

len(x)¶

If x is a dense matrix, returns the product of the number of rows and the number of columns. If x is a sparse matrix, returns the number of nonzero entries.

bool([x])¶

Returns False if x is a zero matrix and True otherwise.

max(x)¶

If x is a dense matrix, returns the maximum element of x. If x is a sparse, returns the maximum nonzero element of x.

min(x)¶

If x is a dense matrix, returns the minimum element of x. If x is a sparse matrix, returns the minimum nonzero element of x.

abs(x)¶

Returns a matrix with the absolute values of the elements of x.

sum(x[, start = 0.0])¶

Returns the sum of start and the elements of x.

Dense and sparse matrices can be used as arguments to the list, tuple, zip, map, and filter functions described in the Python Library Reference. However, one should note that when used with sparse matrix arguments, these functions only consider the nonzero entries. For example, list(A) and tuple(A) construct a list, respectively a tuple, from the elements of A if A is dense, and of the nonzero elements of A if A is sparse.

list(zip(A, B, ...)) returns a list of tuples, with the i-th tuple containing the i-th elements (or nonzero elements) of A, B, ….

>>> from cvxopt import matrix
>>> A = matrix([[-11., -5., -20.], [-6., 0., 7.]])
>>> B = matrix(range(6), (3,2))
>>> list(A)
[-11.0, -5.0, -20.0, -6.0, 0.0, 7.0]
>>> tuple(B)
(0, 1, 2, 3, 4, 5)
>>> list(zip(A, B))
[(-11.0, 0), (-5.0, 1), (-20.0, 2), (-6.0, 3), (0.0, 4), (7.0, 5)]

list(map(f, A)), where f is a function and A is a dense matrix, returns a list constructed by applying f to each element of A. If A is sparse, the function f is applied to each nonzero element of A. Multiple arguments can be provided, for example, as in map(f, A, B), if f is a function with two arguments. In the following example, we return an integer 0-1 matrix with the result of an elementwise comparison.

>>> A = matrix([ [0.5, -0.1, 2.0], [1.5, 0.2, -0.1], [0.3, 1.0, 0.0]])
>>> print(A)
[ 5.00e-01  1.50e+00  3.00e-01]
[-1.00e-01  2.00e-01  1.00e+00]
[ 2.00e+00 -1.00e-01  0.00e+00]
>>> print(matrix(list(map(lambda x: 0 <= x <= 1, A)), A.size))
[ 1  0  1]
[ 0  1  1]
[ 0  0  1]

list(filter(f, A)), where f is a function and A is a matrix, returns a list containing the elements of A (or nonzero elements of A is A is sparse) for which f is true.

>>> A = matrix([[5, -4, 10, -7], [-1, -5, -6, 2], [6, 1, 5, 2],  [-1, 2, -3, -7]])
>>> print(A)
[  5  -1   6  -1]
[ -4  -5   1   2]
[ 10  -6   5  -3]
[ -7   2   2  -7]
>>> list(filter(lambda x: x%2, A))         # list of odd elements in A
[5, -7, -1, -5, 1, 5, -1, -3, -7]
>>> list(filter(lambda x: -2 < x < 3, A))  # list of elements between -2 and 3
[-1, 2, 1, 2, -1, 2]

It is also possible to iterate over matrix elements, as illustrated in the following example.

>>> A = matrix([[5, -3], [9, 11]])
>>> for x in A: print(max(x,0))
...
5
0
9
11
>>> [max(x,0) for x in A]
[5, 0, 9, 11]

The expression x in A returns True if an element of A (or a nonzero element of A if A is sparse) is equal to x and False otherwise.

Other Matrix Functions¶

The following functions can be imported from CVXOPT.

cvxopt.sqrt(x)¶

The elementwise square root of a dense matrix x. The result is returned as a real matrix if x is an integer or real matrix and as a complex matrix if x is a complex matrix. Raises an exception when x is an integer or real matrix with negative elements.

As an example we take the elementwise square root of the sparse matrix

\[\begin{split}A = \left[ \begin{array}{rrrrr} 0 & 2 & 0 & 0 & 3 \\ 2 & 0 & 0 & 0 & 0 \\ 1 & 2 & 0 & 4 & 0 \\ 0 & 0 & 1 & 0 & 0 \end{array} \right]\end{split}\]

>>> from cvxopt import spmatrix, sqrt
>>> A = spmatrix([2,1,2,2,1,3,4], [1,2,0,2,3,0,2], [0,0,1,1,2,3,3])
>>> B = spmatrix(sqrt(A.V), A.I, A.J)
>>> print(B)
[    0      1.41e+00     0      1.73e+00]
[ 1.41e+00     0         0         0    ]
[ 1.00e+00  1.41e+00     0      2.00e+00]
[    0         0      1.00e+00     0    ]

cvxopt.sin(x)¶

The sine function applied elementwise to a dense matrix x. The result is returned as a real matrix if x is an integer or real matrix and as a complex matrix otherwise.

cvxopt.cos(x)¶

The cosine function applied elementwise to a dense matrix x. The result is returned as a real matrix if x is an integer or real matrix and as a complex matrix otherwise.

cvxopt.exp(x)¶

The exponential function applied elementwise to a dense matrix x. The result is returned as a real matrix if x is an integer or real matrix and as a complex matrix otherwise.

cvxopt.log(x)¶

The natural logarithm applied elementwise to a dense matrix x. The result is returned as a real matrix if x is an integer or real matrix and as a complex matrix otherwise. Raises an exception when x is an integer or real matrix with nonpositive elements, or a complex matrix with zero elements.

cvxopt.mul(x0[, x1[, x2 ...]])¶

If the arguments are dense or sparse matrices of the same size, returns the elementwise product of its arguments. The result is a sparse matrix if one or more of its arguments is sparse, and a dense matrix otherwise.

If the arguments include scalars, a scalar product with the scalar is made. (A 1 by 1 dense matrix is treated as a scalar if the dimensions of the other arguments are not all 1 by 1.)

mul can also be called with an iterable (list, tuple, range object, or generator) as its single argument, if the iterable generates a list of dense or sparse matrices or scalars.

>>> from cvxopt import matrix, spmatrix, mul
>>> A = matrix([[1.0, 2.0], [3.0, 4.0]])
>>> B = spmatrix([2.0, 3.0], [0, 1], [0, 1])
>>> print(mul(A, B, -1.0))
[-2.00e+00     0    ]
[    0     -1.20e+01]
>>> print(mul( matrix([k, k+1]) for k in [1,2,3] ))
[  6]
[ 24]

cvxopt.div(x, y)¶

Returns the elementwise division of x by y. x is a dense or sparse matrix, or a scalar (Python number of 1 by 1 dense matrix). y is a dense matrix or a scalar.

cvxopt.max(x0[, x1[, x2 ...]])¶

When called with a single matrix argument, returns the maximum of the elements of the matrix (including the zero entries, if the matrix is sparse).

When called with multiple arguments, the arguments must be matrices of the same size, or scalars, and the elementwise maximum is returned. A 1 by 1 dense matrix is treated as a scalar if the other arguments are not all 1 by 1. If one of the arguments is scalar, and the other arguments are not all 1 by 1, then the scalar argument is interpreted as a dense matrix with all its entries equal to the scalar.

The result is a sparse matrix if all its arguments are sparse matrices. The result is a number if all its arguments are numbers. The result is a dense matrix if at least one of the arguments is a dense matrix.

max can also be called with an iterable (list, tuple, range object, or generator) as its single argument, if the iterable generates a list of dense or sparse matrices or scalars.

>>> from cvxopt import matrix, spmatrix, max
>>> A = spmatrix([2, -3], [0, 1], [0, 1])
>>> print(max(A, -A, 1))
[ 2.00e+00  1.00e+00]
[ 1.00e+00  3.00e+00]

It is important to note the difference between this max and the built-in max, explained in the previous section.

>>> from cvxopt import spmatrix
>>> A = spmatrix([-1.0, -2.0], [0,1], [0,1])
>>> max(A)          # built-in max of a sparse matrix takes maximum over nonzero elements
-1.0
>>> max(A, -1.5)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
NotImplementedError: matrix comparison not implemented
>>> from cvxopt import max
>>> max(A)          # cvxopt.max takes maximum over all the  elements
0.0
>>> print(max(A, -1.5))
[-1.00e+00  0.00e+00]
[ 0.00e+00 -1.50e+00]

cvxopt.min(x0[, x1[, x2 ...]])¶

When called with a single matrix argument, returns the minimum of the elements of the matrix (including the zero entries, if the matrix is sparse).

When called with multiple arguments, the arguments must be matrices of the same size, or scalars, and the elementwise maximum is returned. A 1 by 1 dense matrix is treated as a scalar if the other arguments are not all 1 by 1. If one of the arguments is scalar, and the other arguments are not all 1 by 1, then the scalar argument is interpreted as a dense matrix with all its entries equal to the scalar.

min can also be called with an iterable (list, tuple, range object, or generator) as its single argument, if the iterable generates a list of dense or sparse matrices or scalars.

Randomly Generated Matrices¶

The CVXOPT package provides two functions normal and uniform for generating randomly distributed matrices. The default installation relies on the pseudo-random number generators in the Python standard library random. Alternatively, the random number generators in the GNU Scientific Library (GSL) can be used, if this option is selected during the installation of CVXOPT. The random matrix functions based on GSL are faster than the default functions based on the random module.

cvxopt.normal(nrows[, ncols = 1[, mean = 0.0[, std = 1.0]]])¶

Returns a type 'd' dense matrix of size nrows by ncols with elements chosen from a normal distribution with mean mean and standard deviation std.

cvxopt.uniform(nrows[, ncols = 1[, a = 0.0[, b = 1.0]]])¶

Returns a type 'd' dense matrix of size nrows by ncols matrix with elements uniformly distributed between a and b.

cvxopt.setseed([value])¶

Sets the state of the random number generator. value must be an integer. If value is absent or equal to zero, the value is taken from the system clock. If the Python random number generators are used, this is equivalent to random.seed(value).

cvxopt.getseed()¶

Returns the current state of the random number generator. This function is only available if the GSL random number generators are installed. (The state of the random number generators in the Python random module can be managed via the functions random.getstate and random.setstate.)

Matrix Classes¶

The BLAS exploit several types of matrix structure: symmetric, Hermitian, triangular, and banded. We represent all these matrix classes by dense real or complex matrix objects, with additional arguments that specify the structure.

Vector

A real or complex \(n\)-vector is represented by a matrix of type 'd' or 'z' and length \(n\), with the entries of the vector stored in column-major order.

General matrix

A general real or complex \(m\) by \(n\) matrix is represented by a real or complex matrix of size (\(m\), \(n\)).

Symmetric matrix

A real or complex symmetric matrix of order \(n\) is represented by a real or complex matrix of size (\(n\), \(n\)), and a character argument uplo with two possible values: 'L' and 'U'. If uplo is 'L', the lower triangular part of the symmetric matrix is stored; if uplo is 'U', the upper triangular part is stored. A square matrix X of size (\(n\), \(n\)) can therefore be used to represent the symmetric matrices

\[ \begin{align}\begin{aligned}\begin{split}\left[\begin{array}{ccccc} X[0,0] & X[1,0] & X[2,0] & \cdots & X[n-1,0] \\ X[1,0] & X[1,1] & X[2,1] & \cdots & X[n-1,1] \\ X[2,0] & X[2,1] & X[2,2] & \cdots & X[n-1,2] \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ X[n-1,0] & X[n-1,1] & X[n-1,2] & \cdots & X[n-1,n-1] \end{array}\right] \quad \mbox{(uplo = 'L')},\end{split}\\\begin{split}\left[\begin{array}{ccccc} X[0,0] & X[0,1] & X[0,2] & \cdots & X[0,n-1] \\ X[0,1] & X[1,1] & X[1,2] & \cdots & X[1,n-1] \\ X[0,2] & X[1,2] & X[2,2] & \cdots & X[2,n-1] \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ X[0,n-1] & X[1,n-1] & X[2,n-1] & \cdots & X[n-1,n-1] \end{array}\right] \quad \mbox{(uplo = U')}.\end{split}\end{aligned}\end{align} \]

Complex Hermitian matrix

A complex Hermitian matrix of order \(n\) is represented by a matrix of type 'z' and size (\(n\), \(n\)), and a character argument uplo with the same meaning as for symmetric matrices. A complex matrix X of size (\(n\), \(n\)) can represent the Hermitian matrices

\[ \begin{align}\begin{aligned}\begin{split}\left[\begin{array}{ccccc} \Re X[0,0] & \bar X[1,0] & \bar X[2,0] & \cdots & \bar X[n-1,0] \\ X[1,0] & \Re X[1,1] & \bar X[2,1] & \cdots & \bar X[n-1,1] \\ X[2,0] & X[2,1] & \Re X[2,2] & \cdots & \bar X[n-1,2] \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ X[n-1,0] & X[n-1,1] & X[n-1,2] & \cdots & \Re X[n-1,n-1] \end{array}\right] \quad \mbox{(uplo = 'L')},\end{split}\\\begin{split}\left[\begin{array}{ccccc} \Re X[0,0] & X[0,1] & X[0,2] & \cdots & X[0,n-1] \\ \bar X[0,1] & \Re X[1,1] & X[1,2] & \cdots & X[1,n-1] \\ \bar X[0,2] & \bar X[1,2] & \Re X[2,2] & \cdots & X[2,n-1] \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \bar X[0,n-1] & \bar X[1,n-1] & \bar X[2,n-1] & \cdots & \Re X[n-1,n-1] \end{array}\right] \quad \mbox{(uplo = 'U')}.\end{split}\end{aligned}\end{align} \]

Triangular matrix

A real or complex triangular matrix of order \(n\) is represented by a real or complex matrix of size (\(n\), \(n\)), and two character arguments: an argument uplo with possible values 'L' and 'U' to distinguish between lower and upper triangular matrices, and an argument diag with possible values 'U' and 'N' to distinguish between unit and non-unit triangular matrices. A square matrix X of size (\(n\), \(n\)) can represent the triangular matrices

\[ \begin{align}\begin{aligned}\begin{split}\left[\begin{array}{ccccc} X[0,0] & 0 & 0 & \cdots & 0 \\ X[1,0] & X[1,1] & 0 & \cdots & 0 \\ X[2,0] & X[2,1] & X[2,2] & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ X[n-1,0] & X[n-1,1] & X[n-1,2] & \cdots & X[n-1,n-1] \end{array}\right] \quad \mbox{(uplo = 'L', diag = 'N')},\end{split}\\\begin{split}\left[\begin{array}{ccccc} 1 & 0 & 0 & \cdots & 0 \\ X[1,0] & 1 & 0 & \cdots & 0 \\ X[2,0] & X[2,1] & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ X[n-1,0] & X[n-1,1] & X[n-1,2] & \cdots & 1 \end{array}\right] \quad \mbox{(uplo = 'L', diag = 'U')},\end{split}\\\begin{split}\left[\begin{array}{ccccc} X[0,0] & X[0,1] & X[0,2] & \cdots & X[0,n-1] \\ 0 & X[1,1] & X[1,2] & \cdots & X[1,n-1] \\ 0 & 0 & X[2,2] & \cdots & X[2,n-1] \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & X[n-1,n-1] \end{array}\right] \quad \mbox{(uplo = 'U', diag = 'N')},\end{split}\\\begin{split}\left[\begin{array}{ccccc} 1 & X[0,1] & X[0,2] & \cdots & X[0,n-1] \\ 0 & 1 & X[1,2] & \cdots & X[1,n-1] \\ 0 & 0 & 1 & \cdots & X[2,n-1] \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \end{array}\right] \quad \mbox{(uplo = 'U', diag = 'U')}.\end{split}\end{aligned}\end{align} \]

General band matrix

A general real or complex \(m\) by \(n\) band matrix with \(k_l\) subdiagonals and \(k_u\) superdiagonals is represented by a real or complex matrix X of size (\(k_l + k_u + 1\), \(n\)), and the two integers \(m\) and \(k_l\). The diagonals of the band matrix are stored in the rows of X, starting at the top diagonal, and shifted horizontally so that the entries of column \(k\) of the band matrix are stored in column \(k\) of X. A matrix X of size (\(k_l + k_u + 1\), \(n\)) therefore represents the \(m\) by \(n\) band matrix

\[\begin{split}\left[ \begin{array}{ccccccc} X[k_u,0] & X[k_u-1,1] & X[k_u-2,2] & \cdots & X[0,k_u] & 0 & \cdots \\ X[k_u+1,0] & X[k_u,1] & X[k_u-1,2] & \cdots & X[1,k_u] & X[0,k_u+1] & \cdots \\ X[k_u+2,0] & X[k_u+1,1] & X[k_u,2] & \cdots & X[2,k_u] & X[1,k_u+1] & \cdots \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots \\ X[k_u+k_l,0] & X[k_u+k_l-1,1] & X[k_u+k_l-2,2] & \cdots & & & \\ 0 & X[k_u+k_l,1] & X[k_u+k_l-1,2] & \cdots & & & \\ \vdots & \vdots & \vdots & \ddots & & & \end{array}\right].\end{split}\]

Symmetric band matrix

A real or complex symmetric band matrix of order \(n\) with \(k\) subdiagonals, is represented by a real or complex matrix X of size (\(k+1\), \(n\)), and an argument uplo to indicate whether the subdiagonals (uplo is 'L') or superdiagonals (uplo is 'U') are stored. The \(k+1\) diagonals are stored as rows of X, starting at the top diagonal (i.e., the main diagonal if uplo is 'L', or the \(k\)-th superdiagonal if uplo is 'U') and shifted horizontally so that the entries of the \(k\)-th column of the band matrix are stored in column \(k\) of X. A matrix X of size (\(k+1\), \(n\)) can therefore represent the band matrices

\[ \begin{align}\begin{aligned}\begin{split}\left[ \begin{array}{ccccccc} X[0,0] & X[1,0] & X[2,0] & \cdots & X[k,0] & 0 & \cdots \\ X[1,0] & X[0,1] & X[1,1] & \cdots & X[k-1,1] & X[k,1] & \cdots \\ X[2,0] & X[1,1] & X[0,2] & \cdots & X[k-2,2] & X[k-1,2] & \cdots \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots \\ X[k,0] & X[k-1,1] & X[k-2,2] & \cdots & & & \\ 0 & X[k,1] & X[k-1,2] & \cdots & & & \\ \vdots & \vdots & \vdots & \ddots & & & \end{array}\right] \quad \mbox{(uplo = 'L')},\end{split}\\\begin{split}\left[ \begin{array}{ccccccc} X[k,0] & X[k-1,1] & X[k-2,2] & \cdots & X[0,k] & 0 & \cdots \\ X[k-1,1] & X[k,1] & X[k-1,2] & \cdots & X[1,k] & X[0,k+1] & \cdots \\ X[k-2,2] & X[k-1,2] & X[k,2] & \cdots & X[2,k] & X[1,k+1] & \cdots \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots \\ X[0,k] & X[1,k] & X[2,k] & \cdots & & & \\ 0 & X[0,k+1] & X[1,k+1] & \cdots & & & \\ \vdots & \vdots & \vdots & \ddots & & & \end{array}\right] \quad \mbox{(uplo='U')}.\end{split}\end{aligned}\end{align} \]

Hermitian band matrix

A complex Hermitian band matrix of order \(n\) with \(k\) subdiagonals is represented by a complex matrix of size (\(k+1\), \(n\)) and an argument uplo, with the same meaning as for symmetric band matrices. A matrix X of size (\(k+1\), \(n\)) can represent the band matrices

\[ \begin{align}\begin{aligned}\begin{split}\left[ \begin{array}{ccccccc} \Re X[0,0] & \bar X[1,0] & \bar X[2,0] & \cdots & \bar X[k,0] & 0 & \cdots \\ X[1,0] & \Re X[0,1] & \bar X[1,1] & \cdots & \bar X[k-1,1] & \bar X[k,1] & \cdots \\ X[2,0] & X[1,1] & \Re X[0,2] & \cdots & \bar X[k-2,2] & \bar X[k-1,2] & \cdots \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots \\ X[k,0] & X[k-1,1] & X[k-2,2] & \cdots & & & \\ 0 & X[k,1] & X[k-1,2] & \cdots & & & \\ \vdots & \vdots & \vdots & \ddots & & & \end{array}\right] \quad \mbox{(uplo = 'L')},\end{split}\\\begin{split}\left[ \begin{array}{ccccccc} \Re X[k,0] & X[k-1,1] & X[k-2,2] & \cdots & X[0,k] & 0 & \cdots \\ \bar X[k-1,1] & \Re X[k,1] & X[k-1,2] & \cdots & X[1,k] & X[0,k+1] & \cdots \\ \bar X[k-2,2] & \bar X[k-1,2] & \Re X[k,2] & \cdots & X[2,k] & X[1,k+1] & \cdots \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots \\ \bar X[0,k] & \bar X[1,k] & \bar X[2,k] & \cdots & & & \\ 0 & \bar X[0,k+1] & \bar X[1,k+1] & \cdots & & & \\ \vdots & \vdots & \vdots & \ddots & & & \end{array}\right] \quad \mbox{(uplo='U')}.\end{split}\end{aligned}\end{align} \]

Triangular band matrix

A triangular band matrix of order \(n\) with \(k\) subdiagonals or superdiagonals is represented by a real complex matrix of size (\(k+1\), \(n\)) and two character arguments uplo and diag, with similar conventions as for symmetric band matrices. A matrix X of size (\(k+1\), \(n\)) can represent the band matrices

\[ \begin{align}\begin{aligned}\begin{split}\left[ \begin{array}{cccc} X[0,0] & 0 & 0 & \cdots \\ X[1,0] & X[0,1] & 0 & \cdots \\ X[2,0] & X[1,1] & X[0,2] & \cdots \\ \vdots & \vdots & \vdots & \ddots \\ X[k,0] & X[k-1,1] & X[k-2,2] & \cdots \\ 0 & X[k,1] & X[k-1,1] & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{array}\right] \quad \mbox{(uplo = 'L', diag = 'N')},\end{split}\\\begin{split}\left[ \begin{array}{cccc} 1 & 0 & 0 & \cdots \\ X[1,0] & 1 & 0 & \cdots \\ X[2,0] & X[1,1] & 1 & \cdots \\ \vdots & \vdots & \vdots & \ddots \\ X[k,0] & X[k-1,1] & X[k-2,2] & \cdots \\ 0 & X[k,1] & X[k-1,2] & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{array}\right] \quad \mbox{(uplo = 'L', diag = 'U')},\end{split}\\\begin{split}\left[ \begin{array}{ccccccc} X[k,0] & X[k-1,1] & X[k-2,3] & \cdots & X[0,k] & 0 & \cdots\\ 0 & X[k,1] & X[k-1,2] & \cdots & X[1,k] & X[0,k+1] & \cdots \\ 0 & 0 & X[k,2] & \cdots & X[2,k] & X[1,k+1] & \cdots \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots \end{array}\right] \quad \mbox{(uplo = 'U', diag = 'N')},\end{split}\\\begin{split}\left[ \begin{array}{ccccccc} 1 & X[k-1,1] & X[k-2,3] & \cdots & X[0,k] & 0 & \cdots\\ 0 & 1 & X[k-1,2] & \cdots & X[1,k] & X[0,k+1] & \cdots \\ 0 & 0 & 1 & \cdots & X[2,k] & X[1,k+1] & \cdots \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots \end{array}\right] \quad \mbox{(uplo = 'U', diag = 'U')}.\end{split}\end{aligned}\end{align} \]

When discussing BLAS functions in the following sections we will omit several less important optional arguments that can be used to select submatrices for in-place operations. The complete specification is documented in the docstrings of the source code, and can be viewed with the pydoc help program.

Level 1 BLAS¶

The level 1 functions implement vector operations.

cvxopt.blas.scal(alpha, x)¶

Scales a vector by a constant:

\[x := \alpha x.\]

If x is a real matrix, the scalar argument alpha must be a Python integer or float. If x is complex, alpha can be an integer, float, or complex.

cvxopt.blas.nrm2(x)¶

Euclidean norm of a vector: returns

\[\|x\|_2.\]

cvxopt.blas.asum(x)¶

1-Norm of a vector: returns

\[\|x\|_1 \quad \mbox{($x$ real)}, \qquad \|\Re x\|_1 + \|\Im x\|_1 \quad \mbox{($x$ complex)}.\]

cvxopt.blas.iamax(x)¶

Returns

\[\mathop{\rm argmax}_{k=0,\ldots,n-1} |x_k| \quad \mbox{($x$ real)}, \qquad \mathop{\rm argmax}_{k=0,\ldots,n-1} |\Re x_k| + |\Im x_k| \quad \mbox{($x$ complex)}.\]

If more than one coefficient achieves the maximum, the index of the first \(k\) is returned.

cvxopt.blas.swap(x, y)¶

Interchanges two vectors:

\[x \leftrightarrow y.\]

x and y are matrices of the same type ('d' or 'z').

cvxopt.blas.copy(x, y)¶

Copies a vector to another vector:

\[y := x.\]

x and y are matrices of the same type ('d' or 'z').

cvxopt.blas.axpy(x, y[, alpha = 1.0])¶

Constant times a vector plus a vector:

\[y := \alpha x + y.\]

x and y are matrices of the same type ('d' or 'z'). If x is real, the scalar argument alpha must be a Python integer or float. If x is complex, alpha can be an integer, float, or complex.

cvxopt.blas.dot(x, y)¶

Returns

\[x^Hy.\]

x and y are matrices of the same type ('d' or 'z').

cvxopt.blas.dotu(x, y)¶

Returns

\[x^Ty.\]

x and y are matrices of the same type ('d' or 'z').

Level 2 BLAS¶

The level 2 functions implement matrix-vector products and rank-1 and rank-2 matrix updates. Different types of matrix structure can be exploited using the conventions of the section Matrix Classes.

cvxopt.blas.gemv(A, x, y[, trans = 'N', alpha = 1.0, beta = 0.0])¶

Matrix-vector product with a general matrix:

\[\begin{split}y & := \alpha Ax + \beta y \quad (\mathrm{trans} = \mathrm{'N'}), \\ y & := \alpha A^T x + \beta y \quad (\mathrm{trans} = \mathrm{'T'}), \\ y & := \alpha A^H x + \beta y \quad (\mathrm{trans} = \mathrm{'C'}).\end{split}\]

The arguments A, x, and y must have the same type ('d' or 'z'). Complex values of alpha and beta are only allowed if A is complex.

cvxopt.blas.symv(A, x, y[, uplo = 'L', alpha = 1.0, beta = 0.0])¶

Matrix-vector product with a real symmetric matrix:

\[y := \alpha A x + \beta y,\]

where \(A\) is a real symmetric matrix. The arguments A, x, and y must have type 'd', and alpha and beta must be real.

cvxopt.blas.hemv(A, x, y[, uplo = 'L', alpha = 1.0, beta = 0.0])¶

Matrix-vector product with a real symmetric or complex Hermitian matrix:

\[y := \alpha A x + \beta y,\]

where \(A\) is real symmetric or complex Hermitian. The arguments A, x, y must have the same type ('d' or 'z'). Complex values of alpha and beta are only allowed if A is complex.

cvxopt.blas.trmv(A, x[, uplo = 'L', trans = 'N', diag = 'N'])¶

Matrix-vector product with a triangular matrix:

\[\begin{split}x & := Ax \quad (\mathrm{trans} = \mathrm{'N'}), \\ x & := A^T x \quad (\mathrm{trans} = \mathrm{'T'}), \\ x & := A^H x \quad (\mathrm{trans} = \mathrm{'C'}),\end{split}\]

where \(A\) is square and triangular. The arguments A and x must have the same type ('d' or 'z').

cvxopt.blas.trsv(A, x[, uplo = 'L', trans = 'N', diag = 'N'])¶

Solution of a nonsingular triangular set of linear equations:

\[\begin{split}x & := A^{-1}x \quad (\mathrm{trans} = \mathrm{'N'}), \\ x & := A^{-T}x \quad (\mathrm{trans} = \mathrm{'T'}), \\ x & := A^{-H}x \quad (\mathrm{trans} = \mathrm{'C'}),\end{split}\]

where \(A\) is square and triangular with nonzero diagonal elements. The arguments A and x must have the same type ('d' or 'z').

cvxopt.blas.gbmv(A, m, kl, x, y[, trans = 'N', alpha = 1.0, beta = 0.0])¶

Matrix-vector product with a general band matrix:

\[\begin{split}y & := \alpha Ax + \beta y \quad (\mathrm{trans} = \mathrm{'N'}), \\ y & := \alpha A^T x + \beta y \quad (\mathrm{trans} = \mathrm{'T'}), \\ y & := \alpha A^H x + \beta y \quad (\mathrm{trans} = \mathrm{'C'}),\end{split}\]

where \(A\) is a rectangular band matrix with \(m\) rows and \(k_l\) subdiagonals. The arguments A, x, y must have the same type ('d' or 'z'). Complex values of alpha and beta are only allowed if A is complex.

cvxopt.blas.sbmv(A, x, y[, uplo = 'L', alpha = 1.0, beta = 0.0])¶

Matrix-vector product with a real symmetric band matrix:

\[y := \alpha Ax + \beta y,\]

where \(A\) is a real symmetric band matrix. The arguments A, x, y must have type 'd', and alpha and beta must be real.

cvxopt.blas.hbmv(A, x, y[, uplo = 'L', alpha = 1.0, beta = 0.0])¶

Matrix-vector product with a real symmetric or complex Hermitian band matrix:

\[y := \alpha Ax + \beta y,\]

where \(A\) is a real symmetric or complex Hermitian band matrix. The arguments A, x, y must have the same type ('d' or 'z'). Complex values of alpha and beta are only allowed if A is complex.

cvxopt.blas.tbmv(A, x[, uplo = 'L', trans = 'N', diag = 'N'])¶

Matrix-vector product with a triangular band matrix:

\[\begin{split}x & := Ax \quad (\mathrm{trans} = \mathrm{'N'}), \\ x & := A^T x \quad (\mathrm{trans} = \mathrm{'T'}), \\ x & := A^H x \quad (\mathrm{trans} = \mathrm{'C'}).\end{split}\]

The arguments A and x must have the same type ('d' or 'z').

cvxopt.blas.tbsv(A, x[, uplo = 'L', trans = 'N', diag = 'N'])¶

Solution of a triangular banded set of linear equations:

\[\begin{split}x & := A^{-1}x \quad (\mathrm{trans} = \mathrm{'N'}), \\ x & := A^{-T} x \quad (\mathrm{trans} = \mathrm{'T'}), \\ x & := A^{-H} x \quad (\mathrm{trans} = \mathrm{'T'}),\end{split}\]

where \(A\) is a triangular band matrix of with nonzero diagonal elements. The arguments A and x must have the same type ('d' or 'z').

cvxopt.blas.ger(x, y, A[, alpha = 1.0])¶

General rank-1 update:

\[A := A + \alpha x y^H,\]

where \(A\) is a general matrix. The arguments A, x, and y must have the same type ('d' or 'z'). Complex values of alpha are only allowed if A is complex.

cvxopt.blas.geru(x, y, A[, alpha = 1.0])¶

General rank-1 update:

\[A := A + \alpha x y^T,\]

where \(A\) is a general matrix. The arguments A, x, and y must have the same type ('d' or 'z'). Complex values of alpha are only allowed if A is complex.

cvxopt.blas.syr(x, A[, uplo = 'L', alpha = 1.0])¶

Symmetric rank-1 update:

\[A := A + \alpha xx^T,\]

where \(A\) is a real symmetric matrix. The arguments A and x must have type 'd'. alpha must be a real number.

cvxopt.blas.her(x, A[, uplo = 'L', alpha = 1.0])¶

Hermitian rank-1 update:

\[A := A + \alpha xx^H,\]

where \(A\) is a real symmetric or complex Hermitian matrix. The arguments A and x must have the same type ('d' or 'z'). alpha must be a real number.

cvxopt.blas.syr2(x, y, A[, uplo = 'L', alpha = 1.0])¶

Symmetric rank-2 update:

\[A := A + \alpha (xy^T + yx^T),\]

where \(A\) is a real symmetric matrix. The arguments A, x, and y must have type 'd'. alpha must be real.

cvxopt.blas.her2(x, y, A[, uplo = 'L', alpha = 1.0])¶

Symmetric rank-2 update:

\[A := A + \alpha xy^H + \bar \alpha yx^H,\]

where \(A\) is a a real symmetric or complex Hermitian matrix. The arguments A, x, and y must have the same type ('d' or 'z'). Complex values of alpha are only allowed if A is complex.

As an example, the following code multiplies the tridiagonal matrix

with the vector \(x = (1,-1,2,-2)\).

>>> from cvxopt import matrix
>>> from cvxopt.blas import gbmv
>>> A = matrix([[0., 1., 2.],  [6., -4., -3.],  [3., -1., 0.],  [1., 0., 0.]])
>>> x = matrix([1., -1., 2., -2.])
>>> y = matrix(0., (3,1))
>>> gbmv(A, 3, 1, x, y)
>>> print(y)
[-5.00e+00]
[ 1.20e+01]
[-1.00e+00]

The following example illustrates the use of tbsv.

>>> from cvxopt import matrix
>>> from cvxopt.blas import tbsv
>>> A = matrix([-6., 5., -1., 2.], (1,4))
>>> x = matrix(1.0, (4,1))
>>> tbsv(A, x)  # x := diag(A)^{-1}*x
>>> print(x)
[-1.67e-01]
[ 2.00e-01]
[-1.00e+00]
[ 5.00e-01]

Level 3 BLAS¶

The level 3 BLAS include functions for matrix-matrix multiplication.

cvxopt.blas.gemm(A, B, C[, transA = 'N', transB = 'N', alpha = 1.0, beta = 0.0])¶

Matrix-matrix product of two general matrices:

\[\newcommand{\op}{\mathop{\mathrm{op}}} C := \alpha \op(A) \op(B) + \beta C\]

where

\[\begin{split}\newcommand{\op}{\mathop{\mathrm{op}}} \op(A) = \left\{ \begin{array}{ll} A & \mathrm{transA} = \mathrm{'N'} \\ A^T & \mathrm{transA} = \mathrm{'T'} \\ A^H & \mathrm{transA} = \mathrm{'C'} \end{array} \right. \qquad \op(B) = \left\{ \begin{array}{ll} B & \mathrm{transB} = \mathrm{'N'} \\ B^T & \mathrm{transB} = \mathrm{'T'} \\ B^H & \mathrm{transB} = \mathrm{'C'}. \end{array} \right.\end{split}\]

The arguments A, B, and C must have the same type ('d' or 'z'). Complex values of alpha and beta are only allowed if A is complex.

cvxopt.blas.symm(A, B, C[, side = 'L', uplo = 'L', alpha =1.0, beta = 0.0])¶

Product of a real or complex symmetric matrix \(A\) and a general matrix \(B\):

\[\begin{split}C & := \alpha AB + \beta C \quad (\mathrm{side} = \mathrm{'L'}), \\ C & := \alpha BA + \beta C \quad (\mathrm{side} = \mathrm{'R'}).\end{split}\]

The arguments A, B, and C must have the same type ('d' or 'z'). Complex values of alpha and beta are only allowed if A is complex.

cvxopt.blas.hemm(A, B, C[, side = 'L', uplo = 'L', alpha = 1.0, beta = 0.0])¶

Product of a real symmetric or complex Hermitian matrix \(A\) and a general matrix \(B\):

\[\begin{split}C & := \alpha AB + \beta C \quad (\mathrm{side} = \mathrm{'L'}), \\ C & := \alpha BA + \beta C \quad (\mathrm{side} = \mathrm{'R'}).\end{split}\]

The arguments A, B, and C must have the same type ('d' or 'z'). Complex values of alpha and beta are only allowed if A is complex.

cvxopt.blas.trmm(A, B[, side = 'L', uplo = 'L', transA = 'N', diag = 'N', alpha = 1.0])¶

Product of a triangular matrix \(A\) and a general matrix \(B\):

\[\begin{split}\newcommand{\op}{\mathop{\mathrm{op}}} \begin{split} B & := \alpha\op(A)B \quad (\mathrm{side} = \mathrm{'L'}), \\ B & := \alpha B\op(A) \quad (\mathrm{side} = \mathrm{'R'}) \end{split}\end{split}\]

where

\[\begin{split}\newcommand{\op}{\mathop{\mathrm{op}}} \op(A) = \left\{ \begin{array}{ll} A & \mathrm{transA} = \mathrm{'N'} \\ A^T & \mathrm{transA} = \mathrm{'T'} \\ A^H & \mathrm{transA} = \mathrm{'C'}. \end{array} \right.\end{split}\]

The arguments A and B must have the same type ('d' or 'z'). Complex values of alpha are only allowed if A is complex.

cvxopt.blas.trsm(A, B[, side = 'L', uplo = 'L', transA = 'N', diag = 'N', alpha = 1.0])¶

Solution of a nonsingular triangular system of equations:

\[\begin{split}\newcommand{\op}{\mathop{\mathrm{op}}} \begin{split} B & := \alpha \op(A)^{-1}B \quad (\mathrm{side} = \mathrm{'L'}), \\ B & := \alpha B\op(A)^{-1} \quad (\mathrm{side} = \mathrm{'R'}), \end{split}\end{split}\]

where

\[\begin{split}\newcommand{\op}{\mathop{\mathrm{op}}} \op(A) = \left\{ \begin{array}{ll} A & \mathrm{transA} = \mathrm{'N'} \\ A^T & \mathrm{transA} = \mathrm{'T'} \\ A^H & \mathrm{transA} = \mathrm{'C'}, \end{array} \right.\end{split}\]

\(A\) is triangular and \(B\) is a general matrix. The arguments A and B must have the same type ('d' or 'z'). Complex values of alpha are only allowed if A is complex.

cvxopt.blas.syrk(A, C[, uplo = 'L', trans = 'N', alpha = 1.0, beta = 0.0])¶

Rank-\(k\) update of a real or complex symmetric matrix \(C\):

\[\begin{split}C & := \alpha AA^T + \beta C \quad (\mathrm{trans} = \mathrm{'N'}), \\ C & := \alpha A^TA + \beta C \quad (\mathrm{trans} = \mathrm{'T'}),\end{split}\]

where \(A\) is a general matrix. The arguments A and C must have the same type ('d' or 'z'). Complex values of alpha and beta are only allowed if A is complex.

cvxopt.blas.herk(A, C[, uplo = 'L', trans = 'N', alpha = 1.0, beta = 0.0])¶

Rank-\(k\) update of a real symmetric or complex Hermitian matrix \(C\):

\[\begin{split}C & := \alpha AA^H + \beta C \quad (\mathrm{trans} = \mathrm{'N'}), \\ C & := \alpha A^HA + \beta C \quad (\mathrm{trans} = \mathrm{'C'}),\end{split}\]

where \(A\) is a general matrix. The arguments A and C must have the same type ('d' or 'z'). alpha and beta must be real.

cvxopt.blas.syr2k(A, B, C[, uplo = 'L', trans = 'N', alpha = 1.0, beta = 0.0])¶

Rank-\(2k\) update of a real or complex symmetric matrix \(C\):

\[\begin{split}C & := \alpha (AB^T + BA^T) + \beta C \quad (\mathrm{trans} = \mathrm{'N'}), \\ C & := \alpha (A^TB + B^TA) + \beta C \quad (\mathrm{trans} = \mathrm{'T'}).\end{split}\]

\(A\) and \(B\) are general real or complex matrices. The arguments A, B, and C must have the same type. Complex values of alpha and beta are only allowed if A is complex.

cvxopt.blas.her2k(A, B, C[, uplo = 'L', trans = 'N', alpha = 1.0, beta = 0.0])¶

Rank-\(2k\) update of a real symmetric or complex Hermitian matrix \(C\):

\[\begin{split}C & := \alpha AB^H + \bar \alpha BA^H + \beta C \quad (\mathrm{trans} = \mathrm{'N'}), \\ C & := \alpha A^HB + \bar\alpha B^HA + \beta C \quad (\mathrm{trans} = \mathrm{'C'}),\end{split}\]

where \(A\) and \(B\) are general matrices. The arguments A, B, and C must have the same type ('d' or 'z'). Complex values of alpha are only allowed if A is complex. beta must be real.

General Linear Equations¶

cvxopt.lapack.gesv(A, B[, ipiv = None])¶

Solves

\[A X = B,\]

where \(A\) and \(B\) are real or complex matrices, with \(A\) square and nonsingular.

The arguments A and B must have the same type ('d' or 'z'). On entry, B contains the right-hand side \(B\); on exit it contains the solution \(X\). The optional argument ipiv is an integer matrix of length at least \(n\). If ipiv is provided, then gesv solves the system, replaces A with the triangular factors in an LU factorization, and returns the permutation matrix in ipiv. If ipiv is not specified, then gesv solves the system but does not return the LU factorization and does not modify A.

Raises an ArithmeticError if the matrix is singular.

cvxopt.lapack.getrf(A, ipiv)¶

LU factorization of a general, possibly rectangular, real or complex matrix,

\[A = PLU,\]

where \(A\) is \(m\) by \(n\).

The argument ipiv is an integer matrix of length at least min{\(m\), \(n\)}. On exit, the lower triangular part of A is replaced by \(L\), the upper triangular part by \(U\), and the permutation matrix is returned in ipiv.

Raises an ArithmeticError if the matrix is not full rank.

cvxopt.lapack.getrs(A, ipiv, B[, trans = 'N'])¶

Solves a general set of linear equations

\[\begin{split}AX & = B \quad (\mathrm{trans} = \mathrm{'N'}), \\ A^TX & = B \quad (\mathrm{trans} = \mathrm{'T'}), \\ A^HX & = B \quad (\mathrm{trans} = \mathrm{'C'}),\end{split}\]

given the LU factorization computed by gesv or getrf.

On entry, A and ipiv must contain the factorization as computed by gesv or getrf. On entry, B contains the right-hand side \(B\); on exit it contains the solution \(X\). B must have the same type as A.

cvxopt.lapack.getri(A, ipiv)¶

Computes the inverse of a matrix.

On entry, A and ipiv must contain the factorization as computed by gesv or getrf. On exit, A contains the matrix inverse.

In the following example we compute

for randomly generated problem data, factoring the coefficient matrix once.

>>> from cvxopt import matrix, normal
>>> from cvxopt.lapack import gesv, getrs
>>> n = 10
>>> A = normal(n,n)
>>> b = normal(n)
>>> ipiv = matrix(0, (n,1))
>>> x = +b
>>> gesv(A, x, ipiv)               # x = A^{-1}*b
>>> x2 = +b
>>> getrs(A, ipiv, x2, trans='T')  # x2 = A^{-T}*b
>>> x += x2

Separate functions are provided for equations with band matrices.

cvxopt.lapack.gbsv(A, kl, B[, ipiv = None])¶

Solves

\[A X = B,\]

where \(A\) and \(B\) are real or complex matrices, with \(A\) \(n\) by \(n\) and banded with \(k_l\) subdiagonals.

The arguments A and B must have the same type ('d' or 'z'). On entry, B contains the right-hand side \(B\); on exit it contains the solution \(X\). The optional argument ipiv is an integer matrix of length at least \(n\). If ipiv is provided, then A must have \(2k_l + k_u + 1\) rows. On entry the diagonals of \(A\) are stored in rows \(k_l + 1\) to \(2k_l + k_u + 1\) of A, using the BLAS format for general band matrices (see the section Matrix Classes). On exit, the factorization is returned in A and ipiv. If ipiv is not provided, then A must have \(k_l + k_u + 1\) rows. On entry the diagonals of \(A\) are stored in the rows of A, following the standard BLAS format for general band matrices. In this case, gbsv does not modify A and does not return the factorization.

Raises an ArithmeticError if the matrix is singular.

cvxopt.lapack.gbtrf(A, m, kl, ipiv)¶

LU factorization of a general \(m\) by \(n\) real or complex band matrix with \(k_l\) subdiagonals.

The matrix is stored using the BLAS format for general band matrices (see the section Matrix Classes), by providing the diagonals (stored as rows of a \(k_u + k_l + 1\) by \(n\) matrix A), the number of rows \(m\), and the number of subdiagonals \(k_l\). The argument ipiv is an integer matrix of length at least min{\(m\), \(n\)}. On exit, A and ipiv contain the details of the factorization.

Raises an ArithmeticError if the matrix is not full rank.

cvxopt.lapack.gbtrs({A, kl, ipiv, B[, trans = 'N'])¶

Solves a set of linear equations

\[\begin{split}AX & = B \quad (\mathrm{trans} = \mathrm{'N'}), \\ A^TX & = B \quad (\mathrm{trans} = \mathrm{'T'}), \\ A^HX & = B \quad (\mathrm{trans} = \mathrm{'C'}),\end{split}\]

with \(A\) a general band matrix with \(k_l\) subdiagonals, given the LU factorization computed by gbsv or gbtrf.

On entry, A and ipiv must contain the factorization as computed by gbsv or gbtrf. On entry, B contains the right-hand side \(B\); on exit it contains the solution \(X\). B must have the same type as A.

As an example, we solve a linear equation with

>>> from cvxopt import matrix
>>> from cvxopt.lapack import gbsv, gbtrf, gbtrs
>>> n, kl, ku = 4, 2, 1
>>> A = matrix([[0., 1., 3., 6.], [2., 4., 7., 10.], [5., 8., 11., 0.], [9., 12., 0., 0.]])
>>> x = matrix(1.0, (n,1))
>>> gbsv(A, kl, x)
>>> print(x)
[ 7.14e-02]
[ 4.64e-01]
[-2.14e-01]
[-1.07e-01]

The code below illustrates how one can reuse the factorization returned by gbsv.

>>> Ac = matrix(0.0, (2*kl+ku+1,n))
>>> Ac[kl:,:] = A
>>> ipiv = matrix(0, (n,1))
>>> x = matrix(1.0, (n,1))
>>> gbsv(Ac, kl, x, ipiv)                 # solves A*x = 1
>>> print(x)
[ 7.14e-02]
[ 4.64e-01]
[-2.14e-01]
[-1.07e-01]
>>> x = matrix(1.0, (n,1))
>>> gbtrs(Ac, kl, ipiv, x, trans='T')     # solve A^T*x = 1
>>> print(x)
[ 7.14e-02]
[ 2.38e-02]
[ 1.43e-01]
[-2.38e-02]

An alternative method uses gbtrf for the factorization.

>>> Ac[kl:,:] = A
>>> gbtrf(Ac, n, kl, ipiv)
>>> x = matrix(1.0, (n,1))
>>> gbtrs(Ac, kl, ipiv, x)                # solve A^T*x = 1
>>> print(x)
[ 7.14e-02]
[ 4.64e-01]
[-2.14e-01]
[-1.07e-01]
>>> x = matrix(1.0, (n,1))
>>> gbtrs(Ac, kl, ipiv, x, trans='T')     # solve A^T*x = 1
>>> print(x)
[ 7.14e-02]
[ 2.38e-02]
[ 1.43e-01]
[-2.38e-02]

The following functions can be used for tridiagonal matrices. They use a simpler matrix format, with the diagonals stored in three separate vectors.

cvxopt.lapack.gtsv(dl, d, du, B))¶

Solves

\[A X = B,\]

where \(A\) is an \(n\) by \(n\) tridiagonal matrix.

The subdiagonal of \(A\) is stored as a matrix dl of length \(n-1\), the diagonal is stored as a matrix d of length \(n\), and the superdiagonal is stored as a matrix du of length \(n-1\). The four arguments must have the same type ('d' or 'z'). On exit dl, d, du are overwritten with the details of the LU factorization of \(A\). On entry, B contains the right-hand side \(B\); on exit it contains the solution \(X\).

Raises an ArithmeticError if the matrix is singular.

cvxopt.lapack.gttrf(dl, d, du, du2, ipiv)¶

LU factorization of an \(n\) by \(n\) tridiagonal matrix.

The subdiagonal of \(A\) is stored as a matrix dl of length \(n-1\), the diagonal is stored as a matrix d of length \(n\), and the superdiagonal is stored as a matrix du of length \(n-1\). dl, d and du must have the same type. du2 is a matrix of length \(n-2\), and of the same type as dl. ipiv is an 'i' matrix of length \(n\). On exit, the five arguments contain the details of the factorization.

Raises an ArithmeticError if the matrix is singular.

cvxopt.lapack.gttrs(dl, d, du, du2, ipiv, B[, trans = 'N'])¶

Solves a set of linear equations

\[\begin{split}AX & = B \quad (\mathrm{trans} = \mathrm{'N'}), \\ A^TX & = B \quad (\mathrm{trans} = \mathrm{'T'}), \\ A^HX & = B \quad (\mathrm{trans} = \mathrm{'C'}),\end{split}\]

where \(A\) is an \(n\) by \(n\) tridiagonal matrix.

The arguments dl, d, du, du2, and ipiv contain the details of the LU factorization as returned by gttrf. On entry, B contains the right-hand side \(B\); on exit it contains the solution \(X\). B must have the same type as the other arguments.

Positive Definite Linear Equations¶

cvxopt.lapack.posv(A, B[, uplo = 'L'])¶

Solves

\[A X = B,\]

where \(A\) is a real symmetric or complex Hermitian positive definite matrix.

On exit, B is replaced by the solution, and A is overwritten with the Cholesky factor. The matrices A and B must have the same type ('d' or 'z').

Raises an ArithmeticError if the matrix is not positive definite.

cvxopt.lapack.potrf(A[, uplo = 'L'])¶

Cholesky factorization

\[A = LL^T \qquad \mbox{or} \qquad A = LL^H\]

of a positive definite real symmetric or complex Hermitian matrix \(A\).

On exit, the lower triangular part of A (if uplo is 'L') or the upper triangular part (if uplo is 'U') is overwritten with the Cholesky factor or its (conjugate) transpose.

Raises an ArithmeticError if the matrix is not positive definite.

cvxopt.lapack.potrs(A, B[, uplo = 'L'])¶

Solves a set of linear equations

\[AX = B\]

with a positive definite real symmetric or complex Hermitian matrix, given the Cholesky factorization computed by posv or potrf.

On entry, A contains the triangular factor, as computed by posv or potrf. On exit, B is replaced by the solution. B must have the same type as A.

cvxopt.lapack.potri(A[, uplo = 'L'])¶

Computes the inverse of a positive definite matrix.

On entry, A contains the Cholesky factorization computed by potrf or posv. On exit, it contains the matrix inverse.

As an example, we use posv to solve the linear system

by block-elimination. We first pick a random problem.

>>> from cvxopt import matrix, div, normal, uniform
>>> from cvxopt.blas import syrk, gemv
>>> from cvxopt.lapack import posv
>>> m, n = 100, 50
>>> A = normal(m,n)
>>> b1, b2 = normal(m), normal(n)
>>> d = uniform(m)

We then solve the equations

>>> Asc = div(A, d[:, n*[0]])                # Asc := diag(d)^{-1}*A
>>> B = matrix(0.0, (n,n))
>>> syrk(Asc, B, trans='T')                  # B := Asc^T * Asc = A^T * diag(d)^{-2} * A
>>> x1 = div(b1, d)                          # x1 := diag(d)^{-1}*b1
>>> x2 = +b2
>>> gemv(Asc, x1, x2, trans='T', beta=1.0)   # x2 := x2 + Asc^T*x1 = b2 + A^T*diag(d)^{-2}*b1
>>> posv(B, x2)                              # x2 := B^{-1}*x2 = B^{-1}*(b2 + A^T*diag(d)^{-2}*b1)
>>> gemv(Asc, x2, x1, beta=-1.0)             # x1 := Asc*x2 - x1 = diag(d)^{-1} * (A*x2 - b1)
>>> x1 = div(x1, d)                          # x1 := diag(d)^{-1}*x1 = diag(d)^{-2} * (A*x2 - b1)

There are separate routines for equations with positive definite band matrices.

cvxopt.lapack.pbsv(A, B[, uplo='L'])¶

Solves

\[AX = B\]

where \(A\) is a real symmetric or complex Hermitian positive definite band matrix.

On entry, the diagonals of \(A\) are stored in A, using the BLAS format for symmetric or Hermitian band matrices (see section Matrix Classes). On exit, B is replaced by the solution, and A is overwritten with the Cholesky factor (in the BLAS format for triangular band matrices). The matrices A and B must have the same type ('d' or 'z').

Raises an ArithmeticError if the matrix is not positive definite.

cvxopt.lapack.pbtrf(A[, uplo = 'L'])¶

Cholesky factorization

\[A = LL^T \qquad \mbox{or} \qquad A = LL^H\]

of a positive definite real symmetric or complex Hermitian band matrix \(A\).

On entry, the diagonals of \(A\) are stored in A, using the BLAS format for symmetric or Hermitian band matrices. On exit, A contains the Cholesky factor, in the BLAS format for triangular band matrices.

Raises an ArithmeticError if the matrix is not positive definite.

cvxopt.lapack.pbtrs(A, B[, uplo = 'L'])¶

Solves a set of linear equations

\[AX=B\]

with a positive definite real symmetric or complex Hermitian band matrix, given the Cholesky factorization computed by pbsv or pbtrf.

On entry, A contains the triangular factor, as computed by pbsv or pbtrf. On exit, B is replaced by the solution. B must have the same type as A.

The following functions are useful for tridiagonal systems.

cvxopt.lapack.ptsv(d, e, B)¶

Solves

\[A X = B,\]

where \(A\) is an \(n\) by \(n\) positive definite real symmetric or complex Hermitian tridiagonal matrix.

The diagonal of \(A\) is stored as a 'd' matrix d of length \(n\) and its subdiagonal as a 'd' or 'z' matrix e of length \(n-1\). The arguments e and B must have the same type. On exit d contains the diagonal elements of \(D\) in the LDL^T or LDL^H factorization of \(A\), and e contains the subdiagonal elements of the unit lower bidiagonal matrix \(L\). B is overwritten with the solution \(X\). Raises an ArithmeticError if the matrix is singular.

cvxopt.lapack.pttrf(d, e)¶

LDL^T or LDL^H factorization of an \(n\) by \(n\) positive definite real symmetric or complex Hermitian tridiagonal matrix \(A\).

On entry, the argument d is a 'd' matrix with the diagonal elements of \(A\). The argument e is 'd' or 'z' matrix containing the subdiagonal of \(A\). On exit d contains the diagonal elements of \(D\), and e contains the subdiagonal elements of the unit lower bidiagonal matrix \(L\).

Raises an ArithmeticError if the matrix is singular.

cvxopt.lapack.pttrs(d, e, B[, uplo = 'L'])¶

Solves a set of linear equations

\[AX = B\]

where \(A\) is an \(n\) by \(n\) positive definite real symmetric or complex Hermitian tridiagonal matrix, given its LDL^T or LDL^H factorization.

The argument d is the diagonal of the diagonal matrix \(D\). The argument uplo only matters for complex matrices. If uplo is 'L', then on exit e contains the subdiagonal elements of the unit bidiagonal matrix \(L\). If uplo is 'U', then e contains the complex conjugates of the elements of the unit bidiagonal matrix \(L\). On exit, B is overwritten with the solution \(X\). B must have the same type as e.

Symmetric and Hermitian Linear Equations¶

cvxopt.lapack.sysv(A, B[, ipiv = None, uplo = 'L'])¶

Solves

\[AX = B\]

where \(A\) is a real or complex symmetric matrix of order \(n\).

On exit, B is replaced by the solution. The matrices A and B must have the same type ('d' or 'z'). The optional argument ipiv is an integer matrix of length at least equal to \(n\). If ipiv is provided, sysv solves the system and returns the factorization in A and ipiv. If ipiv is not specified, sysv solves the system but does not return the factorization and does not modify A.

Raises an ArithmeticError if the matrix is singular.

cvxopt.lapack.sytrf(A, ipiv[, uplo = 'L'])¶

LDL^T factorization

\[PAP^T = LDL^T\]

of a real or complex symmetric matrix \(A\) of order \(n\).

ipiv is an 'i' matrix of length at least \(n\). On exit, A and ipiv contain the factorization.

Raises an ArithmeticError if the matrix is singular.

cvxopt.lapack.sytrs(A, ipiv, B[, uplo = 'L'])¶

Solves

\[A X = B\]

given the LDL^T factorization computed by sytrf or sysv. B must have the same type as A.

cvxopt.lapack.sytri(A, ipiv[, uplo = 'L'])¶

Computes the inverse of a real or complex symmetric matrix.

On entry, A and ipiv contain the LDL^T factorization computed by sytrf or sysv. On exit, A contains the inverse.

cvxopt.lapack.hesv(A, B[, ipiv = None, uplo = 'L'])¶

Solves

\[A X = B\]

where \(A\) is a real symmetric or complex Hermitian of order \(n\).

On exit, B is replaced by the solution. The matrices A and B must have the same type ('d' or 'z'). The optional argument ipiv is an integer matrix of length at least \(n\). If ipiv is provided, then hesv solves the system and returns the factorization in A and ipiv. If ipiv is not specified, then hesv solves the system but does not return the factorization and does not modify A.

Raises an ArithmeticError if the matrix is singular.

cvxopt.lapack.hetrf(A, ipiv[, uplo = 'L'])¶

LDL^H factorization

\[PAP^T = LDL^H\]

of a real symmetric or complex Hermitian matrix of order \(n\). ipiv is an 'i' matrix of length at least \(n\). On exit, A and ipiv contain the factorization.

Raises an ArithmeticError if the matrix is singular.

cvxopt.lapack.hetrs(A, ipiv, B[, uplo = 'L'])¶

Solves

\[A X = B\]

given the LDL^H factorization computed by hetrf or hesv.

cvxopt.lapack.hetri(A, ipiv[, uplo = 'L'])¶

Computes the inverse of a real symmetric or complex Hermitian matrix.

On entry, A and ipiv contain the LDL^H factorization computed by hetrf or hesv. On exit, A contains the inverse.

As an example we solve the KKT system (1).

>>> from cvxopt.lapack import sysv
>>> K = matrix(0.0, (m+n,m+n))
>>> K[: (m+n)*m : m+n+1] = -d**2
>>> K[:m, m:] = A
>>> x = matrix(0.0, (m+n,1))
>>> x[:m], x[m:] = b1, b2
>>> sysv(K, x, uplo='U')

Triangular Linear Equations¶

cvxopt.lapack.trtrs(A, B[, uplo = 'L', trans = 'N', diag = 'N'])¶

Solves a triangular set of equations

\[\begin{split}AX & = B \quad (\mathrm{trans} = \mathrm{'N'}), \\ A^TX & = B \quad (\mathrm{trans} = \mathrm{'T'}), \\ A^HX & = B \quad (\mathrm{trans} = \mathrm{'C'}),\end{split}\]

where \(A\) is real or complex and triangular of order \(n\), and \(B\) is a matrix with \(n\) rows.

A and B are matrices with the same type ('d' or 'z'). trtrs is similar to blas.trsm, except that it raises an ArithmeticError if a diagonal element of A is zero (whereas blas.trsm returns inf values).

cvxopt.lapack.trtri(A[, uplo = 'L', diag = 'N'])¶

Computes the inverse of a real or complex triangular matrix \(A\). On exit, A contains the inverse.

cvxopt.lapack.tbtrs(A, B[, uplo = 'L', trans = 'T', diag = 'N'])¶

Solves a triangular set of equations

\[\begin{split}AX & = B \quad (\mathrm{trans} = \mathrm{'N'}), \\ A^TX & = B \quad (\mathrm{trans} = \mathrm{'T'}), \\ A^HX & = B \quad (\mathrm{trans} = \mathrm{'C'}),\end{split}\]

where \(A\) is real or complex triangular band matrix of order \(n\), and \(B\) is a matrix with \(n\) rows.

The diagonals of \(A\) are stored in A using the BLAS conventions for triangular band matrices. A and B are matrices with the same type ('d' or 'z'). On exit, B is replaced by the solution \(X\).

Least-Squares and Least-Norm Problems¶

cvxopt.lapack.gels(A, B[, trans = 'N'])¶

Solves least-squares and least-norm problems with a full rank \(m\) by \(n\) matrix \(A\).

1. trans is 'N'. If \(m\) is greater than or equal to \(n\), gels solves the least-squares problem

   \[\begin{array}{ll} \mbox{minimize} & \|AX-B\|_F. \end{array}\]

   If \(m\) is less than or equal to \(n\), gels solves the least-norm problem

   \[\begin{split}\begin{array}{ll} \mbox{minimize} & \|X\|_F \\ \mbox{subject to} & AX = B. \end{array}\end{split}\]

2. trans is 'T' or 'C' and A and B are real. If \(m\) is greater than or equal to \(n\), gels solves the least-norm problem

   \[\begin{split}\begin{array}{ll} \mbox{minimize} & \|X\|_F \\ \mbox{subject to} & A^TX=B. \end{array}\end{split}\]

   If \(m\) is less than or equal to \(n\), gels solves the least-squares problem

   \[\begin{array}{ll} \mbox{minimize} & \|A^TX-B\|_F. \end{array}\]

3. trans is 'C' and A and B are complex. If \(m\) is greater than or equal to \(n\), gels solves the least-norm problem

   \[\begin{split}\begin{array}{ll} \mbox{minimize} & \|X\|_F \\ \mbox{subject to} & A^HX=B. \end{array}\end{split}\]

   If \(m\) is less than or equal to \(n\), gels solves the least-squares problem

   \[\begin{array}{ll} \mbox{minimize} & \|A^HX-B\|_F. \end{array}\]

A and B must have the same typecode ('d' or 'z'). trans = 'T' is not allowed if A is complex. On exit, the solution \(X\) is stored as the leading submatrix of B. The matrix A is overwritten with details of the QR or the LQ factorization of \(A\).

Note that gels does not check whether \(A\) is full rank.

The following functions compute QR and LQ factorizations.

cvxopt.lapack.geqrf(A, tau)¶

QR factorization of a real or complex matrix A:

\[A = Q R.\]

If \(A\) is \(m\) by \(n\), then \(Q\) is \(m\) by \(m\) and orthogonal/unitary, and \(R\) is \(m\) by \(n\) and upper triangular (if \(m\) is greater than or equal to \(n\)), or upper trapezoidal (if \(m\) is less than or equal to \(n\)).

tau is a matrix of the same type as A and of length min{\(m\), \(n\)}. On exit, \(R\) is stored in the upper triangular/trapezoidal part of A. The matrix \(Q\) is stored as a product of min{\(m\), \(n\)} elementary reflectors in the first min{\(m\), \(n\)} columns of A and in tau.

cvxopt.lapack.gelqf(A, tau)¶

LQ factorization of a real or complex matrix A:

\[A = L Q.\]

If \(A\) is \(m\) by \(n\), then \(Q\) is \(n\) by \(n\) and orthogonal/unitary, and \(L\) is \(m\) by \(n\) and lower triangular (if \(m\) is less than or equal to \(n\)), or lower trapezoidal (if \(m\) is greater than or equal to \(n\)).

tau is a matrix of the same type as A and of length min{\(m\), \(n\)}. On exit, \(L\) is stored in the lower triangular/trapezoidal part of A. The matrix \(Q\) is stored as a product of min{\(m\), \(n\)} elementary reflectors in the first min{\(m\), \(n\)} rows of A and in tau.

cvxopt.lapack.geqp3(A, jpvt, tau)¶

QR factorization with column pivoting of a real or complex matrix \(A\):

\[A P = Q R.\]

If \(A\) is \(m\) by \(n\), then \(Q\) is \(m\) by \(m\) and orthogonal/unitary, and \(R\) is \(m\) by \(n\) and upper triangular (if \(m\) is greater than or equal to \(n\)), or upper trapezoidal (if \(m\) is less than or equal to \(n\)).

tau is a matrix of the same type as A and of length min{\(m\), \(n\)}. jpvt is an integer matrix of length \(n\). On entry, if jpvt[k] is nonzero, then column \(k\) of \(A\) is permuted to the front of \(AP\). Otherwise, column \(k\) is a free column.

On exit, jpvt contains the permutation \(P\): the operation \(AP\) is equivalent to A[:, jpvt-1]. \(R\) is stored in the upper triangular/trapezoidal part of A. The matrix \(Q\) is stored as a product of min{\(m\), \(n\)} elementary reflectors in the first min{\(m\),:math:n} columns of A and in tau.

In most applications, the matrix \(Q\) is not needed explicitly, and it is sufficient to be able to make products with \(Q\) or its transpose. The functions unmqr and ormqr multiply a matrix with the orthogonal matrix computed by geqrf.

cvxopt.lapack.unmqr(A, tau, C[, side = 'L', trans = 'N'])¶

Product with a real orthogonal or complex unitary matrix:

\[\begin{split}\newcommand{\op}{\mathop{\mathrm{op}}} \begin{split} C & := \op(Q)C \quad (\mathrm{side} = \mathrm{'L'}), \\ C & := C\op(Q) \quad (\mathrm{side} = \mathrm{'R'}), \\ \end{split}\end{split}\]

where

\[\begin{split}\newcommand{\op}{\mathop{\mathrm{op}}} \op(Q) = \left\{ \begin{array}{ll} Q & \mathrm{trans} = \mathrm{'N'} \\ Q^T & \mathrm{trans} = \mathrm{'T'} \\ Q^H & \mathrm{trans} = \mathrm{'C'}. \end{array}\right.\end{split}\]

If A is \(m\) by \(n\), then \(Q\) is square of order \(m\) and orthogonal or unitary. \(Q\) is stored in the first min{\(m\), \(n\)} columns of A and in tau as a product of min{\(m\), \(n\)} elementary reflectors, as computed by geqrf. The matrices A, tau, and C must have the same type. trans = 'T' is only allowed if the typecode is 'd'.

cvxopt.lapack.ormqr(A, tau, C[, side = 'L', trans = 'N'])¶

Identical to unmqr but works only for real matrices, and the possible values of trans are 'N' and 'T'.

As an example, we solve a least-squares problem by a direct call to gels, and by separate calls to geqrf, ormqr, and trtrs.

>>> from cvxopt import blas, lapack, matrix, normal
>>> m, n = 10, 5
>>> A, b = normal(m,n), normal(m,1)
>>> x1 = +b
>>> lapack.gels(+A, x1)                  # x1[:n] minimizes || A*x - b ||_2
>>> tau = matrix(0.0, (n,1))
>>> lapack.geqrf(A, tau)                 # A = [Q1, Q2] * [R1; 0]
>>> x2 = +b
>>> lapack.ormqr(A, tau, x2, trans='T')  # x2 := [Q1, Q2]' * x2
>>> lapack.trtrs(A[:n,:], x2, uplo='U')  # x2[:n] := R1^{-1} * x2[:n]
>>> blas.nrm2(x1[:n] - x2[:n])
3.0050798580569307e-16

The next two functions make products with the orthogonal matrix computed by gelqf.

cvxopt.lapack.unmlq(A, tau, C[, side = 'L', trans = 'N'])¶

Product with a real orthogonal or complex unitary matrix:

\[\begin{split}\newcommand{\op}{\mathop{\mathrm{op}}} \begin{split} C & := \op(Q)C \quad (\mathrm{side} = \mathrm{'L'}), \\ C & := C\op(Q) \quad (\mathrm{side} = \mathrm{'R'}), \\ \end{split}\end{split}\]

where

\[\begin{split}\newcommand{\op}{\mathop{\mathrm{op}}} \op(Q) = \left\{ \begin{array}{ll} Q & \mathrm{trans} = \mathrm{'N'}, \\ Q^T & \mathrm{trans} = \mathrm{'T'}, \\ Q^H & \mathrm{trans} = \mathrm{'C'}. \end{array}\right.\end{split}\]

If A is \(m\) by \(n\), then \(Q\) is square of order \(n\) and orthogonal or unitary. \(Q\) is stored in the first min{\(m\), \(n\)} rows of A and in tau as a product of min{\(m\), \(n\)} elementary reflectors, as computed by gelqf. The matrices A, tau, and C must have the same type. trans = 'T' is only allowed if the typecode is 'd'.

cvxopt.lapack.ormlq(A, tau, C[, side = 'L', trans = 'N'])¶

Identical to unmlq but works only for real matrices, and the possible values of trans or 'N' and 'T'.

As an example, we solve a least-norm problem by a direct call to gels, and by separate calls to gelqf, ormlq, and trtrs.

>>> from cvxopt import blas, lapack, matrix, normal
>>> m, n = 5, 10
>>> A, b = normal(m,n), normal(m,1)
>>> x1 = matrix(0.0, (n,1))
>>> x1[:m] = b
>>> lapack.gels(+A, x1)                  # x1 minimizes ||x||_2 subject to A*x = b
>>> tau = matrix(0.0, (m,1))
>>> lapack.gelqf(A, tau)                 # A = [L1, 0] * [Q1; Q2]
>>> x2 = matrix(0.0, (n,1))
>>> x2[:m] = b                           # x2 = [b; 0]
>>> lapack.trtrs(A[:,:m], x2)            # x2[:m] := L1^{-1} * x2[:m]
>>> lapack.ormlq(A, tau, x2, trans='T')  # x2 := [Q1, Q2]' * x2
>>> blas.nrm2(x1 - x2)
0.0