32.18 Dense matrices over the real double field

Module: sage.matrix.matrix_real_double_dense

Dense matrices over the real double field.

Matrix operations use GSl and numpy.

sage: b=Mat(RDF,2,3).basis()
sage: b[0]
[1.0 0.0 0.0]
[0.0 0.0 0.0]

We deal with the case of zero rows or zero columns:

sage: m = MatrixSpace(RDF,0,3)
sage: m.zero_matrix()
[]

TESTS:

sage: a = matrix(RDF,2,range(4), sparse=False)
sage: loads(dumps(a)) == a
True

Author Log:

Class: Matrix_real_double_dense

class Matrix_real_double_dense
Class that implements matrices over the real double field. These are supposed to be fast matrix operations using C doubles. Most operations are implemented using GSl or numpy libraries which will call the underlying BLAS on the system.

sage: m = Matrix(RDF, [[1,2],[3,4]])
sage: m**2
[ 7.0 10.0]
[15.0 22.0]
sage: n= m^(-1); n
[-2.0  1.0]
[ 1.5 -0.5]

To compute eigenvalues the use the functions left eigen_vectors or right_eigenvectors

sage: p,e = m.right_eigenvectors()

the result of eigen is a pair p,e, where p is a list of eigenvalues and the e is a matrix whose columns are the eigenvectors

To solve a linear system Ax = b for A = [[1,2] and b = [5,6] [3,4]]

sage: b = vector(RDF,[5,6])
sage: m.solve_left(b)
(-4.0, 4.5)

See the commands QR,LU,SVD for QR, LU, and singular value decomposition.

Functions: cholesky,$ \,$ determinant,$ \,$ eigenspaces,$ \,$ is_symmetric,$ \,$ left_eigenvectors,$ \,$ log_determinant,$ \,$ LU,$ \,$ LU_valid,$ \,$ numpy,$ \,$ QR,$ \,$ right_eigenvectors,$ \,$ solve_left,$ \,$ solve_left_LU,$ \,$ SVD,$ \,$ transpose

cholesky( )

Return the cholesky factorization of this matrix. The input matrix must be symmetric and positive definite or an exception will be raised.

sage: M = MatrixSpace(RDF,5)
sage: r = matrix(RDF,[[   0.,    0.,    0.,    0.,    1.],[   1.,    1.,    1.,    1.,    1.],[  16.,    8.,    4.,    2.,    1.],[  81.,   27.,    9.,    3.,    1.],[ 256.,   64.,   16.,    4.,    1.]])

sage: m = r*M(1)*r.transpose()
sage: c = m.cholesky()
sage: c*c.transpose()
[ 1.0     1.0     1.0     1.0     1.0]
[ 1.0     5.0    31.0   121.0   341.0]
[ 1.0    31.0   341.0  1555.0  4681.0]
[ 1.0   121.0  1555.0  7381.0 22621.0]
[ 1.0   341.0  4681.0 22621.0 69905.0]

determinant( )

Return the determinant of self.

ALGORITHM: Use GSL (LU decompositon)

sage: m = matrix(RDF,2,range(4)); m.det()
-2.0
sage: m = matrix(RDF,0,[]); m.det()
1.0
sage: m = matrix(RDF, 2, range(6)); m.det()
Traceback (most recent call last):
...
ValueError: self must be square

eigenspaces( )

Return a list of pairs (e, V) where e runs through all complex eigenvalues of this matrix, and V is the corresponding left eigenspace (always a 1-dimensional complex vector space).

sage: m = matrix(RDF, 3, range(9)); m
[0.0 1.0 2.0]
[3.0 4.0 5.0]
[6.0 7.0 8.0]
sage: es = m.eigenspaces()
sage: es # random
[(13.3484692283, Vector space of degree 3 and dimension 1 over Real Double
Field
User basis matrix:
[-0.440242867236 -0.567868371314 -0.695493875393]),
(-1.34846922835, Vector space of degree 3 and dimension 1 over Real Double
Field
User basis matrix:
[-0.897878732262 -0.278434036822  0.341010658618]),
(-9.10854412047e-16, Vector space of degree 3 and dimension 1 over Real
Double Field
User basis matrix:
[ 0.408248290464 -0.816496580928  0.408248290464])]

sage: e, v = es[0]
sage: v = v.basis()[0]
sage: a = v * m
sage: b = e * v
sage: diff = a.change_ring(CDF) - b
sage: abs(abs(diff)) < 1e-10
True
sage: diff # random -- very small numbers
(-2.6645352591e-15, -7.1054273576e-15, -3.5527136788e-15)

is_symmetric( )

Return whether this matrix is symmetric, to the given tolerance.

sage: m = matrix(RDF,2,2,range(4)); m
[0.0 1.0]
[2.0 3.0]
sage: m.is_symmetric()
False
sage: m[1,0] = 1.1; m
[0.0 1.0]
[1.1 3.0]
sage: m.is_symmetric()
False

The tolerance inequality is strict:

sage: m.is_symmetric(tol=0.1)
False
sage: m.is_symmetric(tol=0.11)
True

left_eigenvectors( )

Computes the eigenvalues and *left* eigenvectors of this matrix m acting *from the right*. I.e., vectors v such that v*m = lambda*v.

Output:

eigenvalues
- as a list
corresponding eigenvectors
- as an RDF matrix whose rows are the eigenvectors.

sage: m = Matrix(RDF, 3, range(9))
sage: es = m.left_eigenvectors()
sage: es # random-ish platform-dependent output (low order digits)
([13.3484692283, -1.34846922835, -9.10854412047e-16],
[-0.440242867236 -0.567868371314 -0.695493875393]
[-0.897878732262 -0.278434036822  0.341010658618]
[ 0.408248290464 -0.816496580928  0.408248290464])
sage: e, v = es; e = e[0]; v = v[0]
sage: abs(abs(e*v - v*m)) < 1e-10
True

IMPLEMENTATION: Uses numpy.

log_determinant( )

Compute the log of the absolute value of the determinant using GSL (LU decomposition).

NOTE: This is useful if the usual determinant overlows.

sage: m = matrix(RDF,2,2,range(4)); m
[0.0 1.0]
[2.0 3.0]
sage: RDF(log(abs(m.determinant())))
0.69314718056
sage: m.log_determinant()
0.69314718056
sage: m = matrix(RDF,0,0,[]); m
[]
sage: m.log_determinant()
0.0

LU( )

Computes the LU decomposition of a matrix.

For and square matrix A we can find matrices P,L, and U. s.t. P*A = L*U where P is a permutation matrix, L is lower triangular and U is upper triangular.

Output:

P, L, U
- as a tuple

sage: m = matrix(RDF,4,range(16))
sage: P,L,U = m.LU()
sage: P*m
[12.0 13.0 14.0 15.0]
[ 0.0  1.0  2.0  3.0]
[ 8.0  9.0 10.0 11.0]
[ 4.0  5.0  6.0  7.0]        
sage: L*U
[12.0 13.0 14.0 15.0]
[ 0.0  1.0  2.0  3.0]
[ 8.0  9.0 10.0 11.0]
[ 4.0  5.0  6.0  7.0]

LU_valid( )

Returns True if the LU form of this matrix has already been computed.

sage: A= random_matrix(RDF,3) ; A.LU_valid()
False
sage: L,U,P = A.LU()
sage: A.LU_valid()
True

numpy( )

This method returns a copy of the matrix as a numpy array. It uses the numpy C/api so is very fast.

sage: m = matrix(RDF,[[1,2],[3,4]])
sage: n = m.numpy()
sage: import numpy
sage: numpy.linalg.eig(n)
(array([-0.37228132,  5.37228132]), array([[-0.82456484, -0.41597356],
       [ 0.56576746, -0.90937671]]))
sage: m = matrix(RDF, 2, range(6)); m
[0.0 1.0 2.0]
[3.0 4.0 5.0]
sage: m.numpy()
array([[ 0.,  1.,  2.],
       [ 3.,  4.,  5.]])

TESTS:

sage: m = matrix(RDF,0,5,[]); m
[]
sage: m.numpy()
array([], shape=(0, 5), dtype=float64)
sage: m = matrix(RDF,5,0,[]); m
[]
sage: m.numpy()
array([], shape=(5, 0), dtype=float64)

QR( )

Return the Q,R factorization of a real matrix A.

Input:

self
- a real matrix A

Output:
Q, R
- matrices such that A = Q*R such that the columns of Q are orthogonal (i.e., $ Q^t Q = I$ ), and R is upper triangular.

sage: m = matrix(RDF,3,range(12)); m
[ 0.0  1.0  2.0  3.0]
[ 4.0  5.0  6.0  7.0]
[ 8.0  9.0 10.0 11.0]
sage: Q,R = m.QR()
sage: Q*R
[ 0.0  1.0  2.0  3.0]
[ 4.0  5.0  6.0  7.0]
[ 8.0  9.0 10.0 11.0]

Note that the columns of Q will be an orthogonal

sage: Q*Q.transpose()           # slightly random output.
[1.0                   5.55111512313e-17 -1.11022302463e-16]
[ 5.55111512313e-17    1.0               -5.55111512313e-17]
[-1.11022302463e-16    -5.55111512313e-17               1.0]

right_eigenvectors( )

Computes the eigenvalues and *right* eigenvectors of this matrix m acting *from the left*. I.e., vectors v such that m * v = lambda*v, where v is viewed as a column vector.

Output:

eigenvalues
- as a list
corresponding eigenvectors
- as an RDF matrix whose columns are the eigenvectors.

sage: m = Matrix(RDF, 3, range(9))
sage: m.right_eigenvectors()           # random-ish platform-dependent output (low order digits)
([13.3484692283, -1.34846922835, -6.43352295265e-16],
[ 0.164763817282  0.799699663112  0.408248290464]
[ 0.505774475901  0.104205787719 -0.816496580928]
[ 0.846785134519 -0.591288087674  0.408248290464])

IMPLEMENTATION: Uses numpy.

solve_left( )

Solve the equation A*x = b, where

sage: A = matrix(RDF, 3,3, [1,2,5,7.6,2.3,1,1,2,-1]); A
[ 1.0  2.0  5.0]
[ 7.6  2.3  1.0]
[ 1.0  2.0 -1.0]
sage: b = vector(RDF,[1,2,3])
sage: x = A.solve_left(b); x
(-0.113695090439, 1.39018087855, -0.333333333333)
sage: A*x
(1.0, 2.0, 3.0)

TESTS: We test two degenerate cases:

sage: A = matrix(RDF, 0, 3, [])
sage: A.solve_left(vector(RDF,[]))
()
sage: A = matrix(RDF, 3, 0, [])
sage: A.solve_left(vector(RDF,3, [1,2,3]))
(0.0, 0.0, 0.0)

solve_left_LU( )

Solve the equation A*x = b.

Input:

self
- an invertible matrix
b
- a vector

NOTES: This method precomputes and stores the LU decomposition before solving. If many equations of the form Ax=b need to be solved for a singe matrix A, then this method should be used instead of solve. The first time this method is called it will compute the LU decomposition. If the matrix has not changed then subsequent calls will be very fast as the precomputed LU decomposition will be reused.

sage: A = matrix(RDF, 3,3, [1,2,5,7.6,2.3,1,1,2,-1]); A
[ 1.0  2.0  5.0]
[ 7.6  2.3  1.0]
[ 1.0  2.0 -1.0]
sage: b = vector(RDF,[1,2,3])
sage: x = A.solve_left(b); x
(-0.113695090439, 1.39018087855, -0.333333333333)
sage: A*x
(1.0, 2.0, 3.0)

TESTS: We test two degenerate cases:

sage: A = matrix(RDF, 0, 3, [])
sage: A.solve_left_LU(vector(RDF,[]))
(0.0, 0.0, 0.0)
sage: A = matrix(RDF, 3, 0, [])
sage: A.solve_left_LU(vector(RDF,3, [1,2,3]))
()

SVD( )

Return the singular value decomposition of this matrix.

Input:

A
- a matrix
algorithm
- 'numpy' or 'gsl'
Output:
U, S, V
- matrices such that $ A = U S V^t$ , where U and V are orthogonal and S is diagonal.

sage: m = matrix(RDF,4,range(16))
sage: U,S,V = m.SVD(algorithm='gsl')
sage: U*S*V.transpose()    # slightly random output (due to computer architecture)
[3.45569519412e-16               1.0               2.0               3.0]
[4.0               5.0               6.0               7.0]
[8.0               9.0              10.0              11.0]
[12.0              13.0              14.0              15.0]

A non-square example:

sage: m = matrix(RDF, 2, range(6)); m
[0.0 1.0 2.0]
[3.0 4.0 5.0]
sage: U, S, V = m.SVD(algorithm='numpy')
sage: U*S*V.transpose()           # random low bits
[7.62194127257e-17               1.0               2.0]
[              3.0               4.0               5.0]            
sage: U, S, V = m.SVD()
sage: U
[-0.274721127897 -0.961523947641]
[-0.961523947641  0.274721127897]
sage: S
[7.34846922835           0.0]
[          0.0           1.0]
sage: V
[-0.392540507864  0.824163383692]
[-0.560772154092  0.137360563949]
[ -0.72900380032 -0.549442255795]
sage: U*S*V.transpose()           # random low bits
[7.62194127257e-17               1.0               2.0]
[              3.0               4.0               5.0]
sage: m = matrix(RDF,3,2,range(6)); m
[0.0 1.0]
[2.0 3.0]
[4.0 5.0]
sage: U,S,V = m.SVD(algorithm='numpy')
sage: U*S*V.transpose()   # random low order bits
[-8.13151629364e-19                1.0]
[               2.0                3.0]
[               4.0                5.0]  
sage: U,S,V = m.SVD()
sage: U*S*V.transpose()   # random low order bits
[-8.13151629364e-19                1.0]
[               2.0                3.0]
[               4.0                5.0]

TESTS:

sage: m = matrix(RDF, 3, 0, []); m
[]
sage: m.SVD()
([], [], [])
sage: m = matrix(RDF, 0, 3, []); m
[]        
sage: m.SVD()
([], [], [])

transpose( )

Return the transpose of this matrix.

sage: m = matrix(RDF,2,3,range(6)); m
[0.0 1.0 2.0]
[3.0 4.0 5.0]
sage: m.transpose()
[0.0 3.0]
[1.0 4.0]
[2.0 5.0]
sage: m = matrix(RDF,0,3); m
[]
sage: m.transpose()
[]
sage: m.transpose().parent()
Full MatrixSpace of 3 by 0 dense matrices over Real Double Field

Special Functions: __eq__,$ \,$ __ge__,$ \,$ __gt__,$ \,$ __init__,$ \,$ __invert__,$ \,$ __le__,$ \,$ __lt__,$ \,$ __ne__,$ \,$ __neg__,$ \,$ _cholesky_gsl,$ \,$ _hadamard_row_bound,$ \,$ _replace_self_with_numpy,$ \,$ _replace_self_with_numpy32,$ \,$ _SVD_gsl,$ \,$ _SVD_numpy

__invert__( )

Invert this matrix.

sage: A = Matrix(RDF, [[10, 0], [0, 100]])
sage: (~A).det()
0.001

sage: A = matrix(RDF,3,[2,3,5,7,8,9,11,13,17]);A
[ 2.0  3.0  5.0]
[ 7.0  8.0  9.0]
[11.0 13.0 17.0]
sage: ~A
[ -2.71428571429            -2.0   1.85714285714]
[  2.85714285714             3.0  -2.42857142857]
[-0.428571428571            -1.0  0.714285714286]

Note that if this matrix is (nearly) singular, finding its inverse will not help much and will give slightly different answers on similar platforms depending on the hardware and tuning options given to ATLAS:

sage: A = Matrix(RDF, [[1, 0], [0, 0]])
sage: A.inverse().det() 
nan
sage: A = matrix(RDF,3,range(1,10));A
[1.0 2.0 3.0]
[4.0 5.0 6.0]
[7.0 8.0 9.0]

sage: A.determinant() < 10e-12
True
sage: ~A              # slightly random
[-4.50359962737e+15  9.00719925474e+15 -4.50359962737e+15]
[ 9.00719925474e+15 -1.80143985095e+16  9.00719925474e+15]
[-4.50359962737e+15  9.00719925474e+15 -4.50359962737e+15]

__neg__( )

Negate this matrix

sage: A = matrix(RDF,3,range(1,10))
sage: -A
[-1.0 -2.0 -3.0]
[-4.0 -5.0 -6.0]
[-7.0 -8.0 -9.0]
sage: B = -A ; (A+B).is_zero()
True

_cholesky_gsl( )

_hadamard_row_bound( )

Return an integer n such that the absolute value of the determinant of this matrix is at most $ 10^n$ .

sage: a = matrix(RDF, 3, [1,2,5,7,-3,4,2,1,123])
sage: a._hadamard_row_bound()
4
sage: a.det()
-2014.0
sage: 10^4
10000

_replace_self_with_numpy( )

sage: import numpy
sage: a = numpy.array([[1,2],[3,4]], 'float64')
sage: m = matrix(RDF,2,2,0)
sage: m._replace_self_with_numpy(a)
sage: m
[1.0 2.0]
[3.0 4.0]

_replace_self_with_numpy32( )

sage: import numpy
sage: a = numpy.array([[1,2],[3,4]], 'float32')
sage: m = matrix(RDF,2,2,0)
sage: m._replace_self_with_numpy32(a)
sage: m
[1.0 2.0]
[3.0 4.0]

_SVD_gsl( )

Return the singular value decomposition of this matrix using GSL. Note that the matrices the dimensions of the matrices that this returns are (m,p), (p,p), and (n, p) where p = min(m,n).

sage: def shape(x): return (x.nrows(), x.ncols())
sage: m = matrix(RDF, 2, 3, range(6))
sage: map(shape, m._SVD_gsl())
[(2, 2), (2, 2), (3, 2)]

_SVD_numpy( )

Return the singular value decomposition of this matrix using GSL. Note that the matrices the dimensions of the matrices that this returns are (m,m), (m,n), and (n, n).

sage: def shape(x): return (x.nrows(), x.ncols())
sage: m = matrix(RDF, 2, 3, range(6))
sage: map(shape, m._SVD_numpy())
[(2, 2), (2, 3), (3, 3)]

See About this document... for information on suggesting changes.