סדנאת פייתון לפיזיקאים
SciPy, אוסף ספריות פייתון לשימוש מדעי
הפקולטה לפיזיקה, הטכניון. חורף 2013
מרצה: רונן אברבנאל
פייתון היא שפת תכנות לשימוש כללי, וככזו, היא לא כוללת כלים נפוצים שהיינו מצפים למצוא בסביבת עבודה כזו. ב- SciPy.org
מרכזים אוסף של ספריות, סטנדרטיות יחסית, המשלימות את החסך הזה.
SciPy (pronounced “Sigh Pie”) is a Python-based ecosystem of open-source software for mathematics, science, and engineering. In particular, these are some of the core packages:
-
NumPy
Base N-dimensional array package -
SciPy library
Fundamental library for scientific computing
-
Matplotlib
Comprehensive 2D Plotting -
IPython
Enhanced Interactive Console -
Sympy
Symbolic mathematics -
pandas
Data structures & analysis
בשום אופן הן לא הספריות היחידות. כמעט לכל אחת והאחת מהספריות הללו קיימת חלופה סבירה. אבל אלו הספריות ה"סטנדרטיות", הנפוצות ביותר.
יוצאת מהכלל כאן הואnumpy, שמגדירה את הסטנדרט בהתייחסות למערכים ומטריצות, וכמעט כל ספריית פייתון מודרנים עם שימוש מדעי/הנדסי/מתמטי כלשהו, משתמשת בה בתור בסיס, או ממשרת את אותו הממשק.
כל הספריות הנ"ל בפיתוח רציף, וגרסאות חדשות, עם תכונות חדשות ותיקוני באגים, משתחררות אחת ל~חצי שנה.
האלמנט הבסיסי ב-numpy הוא המערכים. מערכים ב-numpy דומים לרשימות של פייתון, אבל הטיפול בהם הוא הרבה פעמים יעיל יותר, וניתן לבצע פעולות "מתמטיות" יותר בקלות רבה יותר.
In [6]:
import numpy as np
a = np.array([0,1,2,3,4])
a
Out[6]:
array([0, 1, 2, 3, 4])
In [9]:
L = range(10000)
%timeit [i**2 for i in L]
1000 loops, best of 3: 754 µs per loop
In [10]:
a = np.arange(10000)
%timeit a**2
100000 loops, best of 3: 11.4 µs per loop
בשורות הללו אנחנו יכולים לראות שני דברים: העלאת כל איבר בריבוע (או באופן כללי, ביצוע פעולה חשבונית על כל אחד מהאיברים) פשוטה הרבה יותר לביצוע כשמדובר ב-numpy array והיא מבוצעת (הרבה) יותר מהר.
In [7]:
L = [1,2,3]
L+L
Out[7]:
[1, 2, 3, 1, 2, 3]
In [9]:
a = np.array([1,2,3])
a+a
Out[9]:
array([2, 4, 6])
In [74]:
a = np.random.randint(0,10,(5,5))
print a.shape
a
(5, 5)
Out[74]:
array([[1, 6, 8, 6, 5],
[3, 6, 5, 7, 8],
[4, 1, 2, 8, 7],
[6, 4, 7, 9, 7],
[8, 0, 6, 7, 6]])
In [75]:
a*a
Out[75]:
array([[ 1, 36, 64, 36, 25],
[ 9, 36, 25, 49, 64],
[16, 1, 4, 64, 49],
[36, 16, 49, 81, 49],
[64, 0, 36, 49, 36]])
In [76]:
a.dot(a)
Out[76]:
array([[127, 74, 126, 201, 181],
[147, 87, 161, 219, 195],
[119, 64, 139, 168, 140],
[156, 103, 187, 250, 216],
[122, 82, 161, 201, 167]])
In [77]:
b = np.matrix(a)
b*b
Out[77]:
matrix([[127, 74, 126, 201, 181],
[147, 87, 161, 219, 195],
[119, 64, 139, 168, 140],
[156, 103, 187, 250, 216],
[122, 82, 161, 201, 167]])
המערכים של numpy הם, כאמור, השלד, ועוד נחזור אליו בהמשך.
הספריה Scipy
scipy logo
כאן נדבר על הספריה scipy ולא על האוסף בכללותו. כלומר, לא על כל הספריות שראינו לעיל, אלא פשוט על "מה שמקבלים" כשמבצעים
import scipyל-scipy יש הרבה תת חבילות, כל אחת מהן היא עצמאית וכל אחת מהן שימושית למדי כשלעצמה.
| scipy.cluster | Vector quantization / Kmeans |
| scipy.constants | Physical and mathematical constants |
| scipy.fftpack | Fourier transform |
| scipy.integrate | Integration routines |
| scipy.interpolate | Interpolation |
| scipy.io | Data input and output |
| scipy.linalg | Linear algebra routines |
| scipy.ndimage | n-dimensional image package |
| scipy.odr | Orthogonal distance regression |
| scipy.optimize | Optimization |
| scipy.signal | Signal processing |
| scipy.sparse | Sparse matrices |
| scipy.spatial | Spatial data structures and algorithms |
| scipy.special | Any special mathematical functions |
| scipy.stats | Statistics |
בהרצאה הזו, לא נדון בהן בצורה ריגורוזית, אלא ניתן דוגמאות כשנדרש להן.
בכל מקרה,
In [13]:
import scipy.linalg
help(scipy.linalg)
?scipy.linalg #in IPython.
Help on package scipy.linalg in scipy:
NAME
scipy.linalg
FILE
/usr/lib/python2.7/dist-packages/scipy/linalg/__init__.py
DESCRIPTION
====================================
Linear algebra (:mod:`scipy.linalg`)
====================================
.. currentmodule:: scipy.linalg
Linear algebra functions.
.. seealso::
`numpy.linalg` for more linear algebra functions. Note that
although `scipy.linalg` imports most of them, identically named
functions from `scipy.linalg` may offer more or slightly differing
functionality.
Basics
======
.. autosummary::
:toctree: generated/
inv - Find the inverse of a square matrix
solve - Solve a linear system of equations
solve_banded - Solve a banded linear system
solveh_banded - Solve a Hermitian or symmetric banded system
solve_triangular - Solve a triangular matrix
det - Find the determinant of a square matrix
norm - Matrix and vector norm
lstsq - Solve a linear least-squares problem
pinv - Pseudo-inverse (Moore-Penrose) using lstsq
pinv2 - Pseudo-inverse using svd
kron - Kronecker product of two arrays
tril - Construct a lower-triangular matrix from a given matrix
triu - Construct an upper-triangular matrix from a given matrix
Eigenvalue Problems
===================
.. autosummary::
:toctree: generated/
eig - Find the eigenvalues and eigenvectors of a square matrix
eigvals - Find just the eigenvalues of a square matrix
eigh - Find the e-vals and e-vectors of a Hermitian or symmetric matrix
eigvalsh - Find just the eigenvalues of a Hermitian or symmetric matrix
eig_banded - Find the eigenvalues and eigenvectors of a banded matrix
eigvals_banded - Find just the eigenvalues of a banded matrix
Decompositions
==============
.. autosummary::
:toctree: generated/
lu - LU decomposition of a matrix
lu_factor - LU decomposition returning unordered matrix and pivots
lu_solve - Solve Ax=b using back substitution with output of lu_factor
svd - Singular value decomposition of a matrix
svdvals - Singular values of a matrix
diagsvd - Construct matrix of singular values from output of svd
orth - Construct orthonormal basis for the range of A using svd
cholesky - Cholesky decomposition of a matrix
cholesky_banded - Cholesky decomp. of a sym. or Hermitian banded matrix
cho_factor - Cholesky decomposition for use in solving a linear system
cho_solve - Solve previously factored linear system
cho_solve_banded - Solve previously factored banded linear system
qr - QR decomposition of a matrix
qr_multiply - QR decomposition and multiplication by Q
qz - QZ decomposition of a pair of matrices
schur - Schur decomposition of a matrix
rsf2csf - Real to complex Schur form
hessenberg - Hessenberg form of a matrix
Matrix Functions
================
.. autosummary::
:toctree: generated/
expm - Matrix exponential using Pade approximation
expm2 - Matrix exponential using eigenvalue decomposition
expm3 - Matrix exponential using Taylor-series expansion
logm - Matrix logarithm
cosm - Matrix cosine
sinm - Matrix sine
tanm - Matrix tangent
coshm - Matrix hyperbolic cosine
sinhm - Matrix hyperbolic sine
tanhm - Matrix hyperbolic tangent
signm - Matrix sign
sqrtm - Matrix square root
funm - Evaluating an arbitrary matrix function
Matrix Equation Solvers
=======================
.. autosummary::
:toctree: generated/
solve_sylvester - Solve the Sylvester matrix equation
solve_continuous_are - Solve the continuous-time algebraic Riccati equation
solve_discrete_are - Solve the discrete-time algebraic Riccati equation
solve_discrete_lyapunov - Solve the discrete-time Lyapunov equation
solve_lyapunov - Solve the (continous-time) Lyapunov equation
Special Matrices
================
.. autosummary::
:toctree: generated/
block_diag - Construct a block diagonal matrix from submatrices
circulant - Circulant matrix
companion - Companion matrix
hadamard - Hadamard matrix of order 2**n
hankel - Hankel matrix
hilbert - Hilbert matrix
invhilbert - Inverse Hilbert matrix
leslie - Leslie matrix
pascal - Pascal matrix
toeplitz - Toeplitz matrix
tri - Construct a matrix filled with ones at and below a given diagonal
PACKAGE CONTENTS
_decomp_qz
_flinalg
_solvers
_testutils
atlas_version
basic
blas
calc_lwork
cblas
clapack
decomp
decomp_cholesky
decomp_lu
decomp_qr
decomp_schur
decomp_svd
fblas
flapack
flinalg
interface_gen
lapack
linalg_version
matfuncs
misc
scons_support
setup
setup_atlas_version
setupscons
special_matrices
CLASSES
exceptions.Exception(exceptions.BaseException)
numpy.linalg.linalg.LinAlgError
class LinAlgError(exceptions.Exception)
| Generic Python-exception-derived object raised by linalg functions.
|
| General purpose exception class, derived from Python's exception.Exception
| class, programmatically raised in linalg functions when a Linear
| Algebra-related condition would prevent further correct execution of the
| function.
|
| Parameters
| ----------
| None
|
| Examples
| --------
| >>> from numpy import linalg as LA
| >>> LA.inv(np.zeros((2,2)))
| Traceback (most recent call last):
| File "<stdin>", line 1, in <module>
| File "...linalg.py", line 350,
| in inv return wrap(solve(a, identity(a.shape[0], dtype=a.dtype)))
| File "...linalg.py", line 249,
| in solve
| raise LinAlgError('Singular matrix')
| numpy.linalg.LinAlgError: Singular matrix
|
| Method resolution order:
| LinAlgError
| exceptions.Exception
| exceptions.BaseException
| __builtin__.object
|
| Data descriptors defined here:
|
| __weakref__
| list of weak references to the object (if defined)
|
| ----------------------------------------------------------------------
| Methods inherited from exceptions.Exception:
|
| __init__(...)
| x.__init__(...) initializes x; see help(type(x)) for signature
|
| ----------------------------------------------------------------------
| Data and other attributes inherited from exceptions.Exception:
|
| __new__ = <built-in method __new__ of type object>
| T.__new__(S, ...) -> a new object with type S, a subtype of T
|
| ----------------------------------------------------------------------
| Methods inherited from exceptions.BaseException:
|
| __delattr__(...)
| x.__delattr__('name') <==> del x.name
|
| __getattribute__(...)
| x.__getattribute__('name') <==> x.name
|
| __getitem__(...)
| x.__getitem__(y) <==> x[y]
|
| __getslice__(...)
| x.__getslice__(i, j) <==> x[i:j]
|
| Use of negative indices is not supported.
|
| __reduce__(...)
|
| __repr__(...)
| x.__repr__() <==> repr(x)
|
| __setattr__(...)
| x.__setattr__('name', value) <==> x.name = value
|
| __setstate__(...)
|
| __str__(...)
| x.__str__() <==> str(x)
|
| __unicode__(...)
|
| ----------------------------------------------------------------------
| Data descriptors inherited from exceptions.BaseException:
|
| __dict__
|
| args
|
| message
FUNCTIONS
all_mat(*args)
block_diag(*arrs)
Create a block diagonal matrix from provided arrays.
Given the inputs `A`, `B` and `C`, the output will have these
arrays arranged on the diagonal::
[[A, 0, 0],
[0, B, 0],
[0, 0, C]]
Parameters
----------
A, B, C, ... : array_like, up to 2-D
Input arrays. A 1-D array or array_like sequence of length `n`is
treated as a 2-D array with shape ``(1,n)``.
Returns
-------
D : ndarray
Array with `A`, `B`, `C`, ... on the diagonal. `D` has the
same dtype as `A`.
Notes
-----
If all the input arrays are square, the output is known as a
block diagonal matrix.
Examples
--------
>>> from scipy.linalg import block_diag
>>> A = [[1, 0],
... [0, 1]]
>>> B = [[3, 4, 5],
... [6, 7, 8]]
>>> C = [[7]]
>>> block_diag(A, B, C)
[[1 0 0 0 0 0]
[0 1 0 0 0 0]
[0 0 3 4 5 0]
[0 0 6 7 8 0]
[0 0 0 0 0 7]]
>>> block_diag(1.0, [2, 3], [[4, 5], [6, 7]])
array([[ 1., 0., 0., 0., 0.],
[ 0., 2., 3., 0., 0.],
[ 0., 0., 0., 4., 5.],
[ 0., 0., 0., 6., 7.]])
cho_factor(a, lower=False, overwrite_a=False)
Compute the Cholesky decomposition of a matrix, to use in cho_solve
Returns a matrix containing the Cholesky decomposition,
``A = L L*`` or ``A = U* U`` of a Hermitian positive-definite matrix `a`.
The return value can be directly used as the first parameter to cho_solve.
.. warning::
The returned matrix also contains random data in the entries not
used by the Cholesky decomposition. If you need to zero these
entries, use the function `cholesky` instead.
Parameters
----------
a : array, shape (M, M)
Matrix to be decomposed
lower : boolean
Whether to compute the upper or lower triangular Cholesky factorization
(Default: upper-triangular)
overwrite_a : boolean
Whether to overwrite data in a (may improve performance)
Returns
-------
c : array, shape (M, M)
Matrix whose upper or lower triangle contains the Cholesky factor
of `a`. Other parts of the matrix contain random data.
lower : boolean
Flag indicating whether the factor is in the lower or upper triangle
Raises
------
LinAlgError
Raised if decomposition fails.
See also
--------
cho_solve : Solve a linear set equations using the Cholesky factorization
of a matrix.
cho_solve((c, lower), b, overwrite_b=False)
Solve the linear equations A x = b, given the Cholesky factorization of A.
Parameters
----------
(c, lower) : tuple, (array, bool)
Cholesky factorization of a, as given by cho_factor
b : array
Right-hand side
Returns
-------
x : array
The solution to the system A x = b
See also
--------
cho_factor : Cholesky factorization of a matrix
cho_solve_banded((cb, lower), b, overwrite_b=False)
Solve the linear equations A x = b, given the Cholesky factorization of A.
Parameters
----------
(cb, lower) : tuple, (array, bool)
`cb` is the Cholesky factorization of A, as given by cholesky_banded.
`lower` must be the same value that was given to cholesky_banded.
b : array
Right-hand side
overwrite_b : bool
If True, the function will overwrite the values in `b`.
Returns
-------
x : array
The solution to the system A x = b
See also
--------
cholesky_banded : Cholesky factorization of a banded matrix
Notes
-----
.. versionadded:: 0.8.0
cholesky(a, lower=False, overwrite_a=False)
Compute the Cholesky decomposition of a matrix.
Returns the Cholesky decomposition, :math:`A = L L^*` or
:math:`A = U^* U` of a Hermitian positive-definite matrix A.
Parameters
----------
a : ndarray, shape (M, M)
Matrix to be decomposed
lower : bool
Whether to compute the upper or lower triangular Cholesky
factorization. Default is upper-triangular.
overwrite_a : bool
Whether to overwrite data in `a` (may improve performance).
Returns
-------
c : ndarray, shape (M, M)
Upper- or lower-triangular Cholesky factor of `a`.
Raises
------
LinAlgError : if decomposition fails.
Examples
--------
>>> from scipy import array, linalg, dot
>>> a = array([[1,-2j],[2j,5]])
>>> L = linalg.cholesky(a, lower=True)
>>> L
array([[ 1.+0.j, 0.+0.j],
[ 0.+2.j, 1.+0.j]])
>>> dot(L, L.T.conj())
array([[ 1.+0.j, 0.-2.j],
[ 0.+2.j, 5.+0.j]])
cholesky_banded(ab, overwrite_ab=False, lower=False)
Cholesky decompose a banded Hermitian positive-definite matrix
The matrix a is stored in ab either in lower diagonal or upper
diagonal ordered form:
ab[u + i - j, j] == a[i,j] (if upper form; i <= j)
ab[ i - j, j] == a[i,j] (if lower form; i >= j)
Example of ab (shape of a is (6,6), u=2)::
upper form:
* * a02 a13 a24 a35
* a01 a12 a23 a34 a45
a00 a11 a22 a33 a44 a55
lower form:
a00 a11 a22 a33 a44 a55
a10 a21 a32 a43 a54 *
a20 a31 a42 a53 * *
Parameters
----------
ab : array, shape (u + 1, M)
Banded matrix
overwrite_ab : boolean
Discard data in ab (may enhance performance)
lower : boolean
Is the matrix in the lower form. (Default is upper form)
Returns
-------
c : array, shape (u+1, M)
Cholesky factorization of a, in the same banded format as ab
circulant(c)
Construct a circulant matrix.
Parameters
----------
c : array_like
1-D array, the first column of the matrix.
Returns
-------
A : array, shape (len(c), len(c))
A circulant matrix whose first column is `c`.
See also
--------
toeplitz : Toeplitz matrix
hankel : Hankel matrix
Notes
-----
.. versionadded:: 0.8.0
Examples
--------
>>> from scipy.linalg import circulant
>>> circulant([1, 2, 3])
array([[1, 3, 2],
[2, 1, 3],
[3, 2, 1]])
companion(a)
Create a companion matrix.
Create the companion matrix [1]_ associated with the polynomial whose
coefficients are given in `a`.
Parameters
----------
a : array_like
1-D array of polynomial coefficients. The length of `a` must be
at least two, and ``a[0]`` must not be zero.
Returns
-------
c : ndarray
A square array of shape ``(n-1, n-1)``, where `n` is the length
of `a`. The first row of `c` is ``-a[1:]/a[0]``, and the first
sub-diagonal is all ones. The data-type of the array is the same
as the data-type of ``1.0*a[0]``.
Raises
------
ValueError
If any of the following are true: a) ``a.ndim != 1``;
b) ``a.size < 2``; c) ``a[0] == 0``.
Notes
-----
.. versionadded:: 0.8.0
References
----------
.. [1] R. A. Horn & C. R. Johnson, *Matrix Analysis*. Cambridge, UK:
Cambridge University Press, 1999, pp. 146-7.
Examples
--------
>>> from scipy.linalg import companion
>>> companion([1, -10, 31, -30])
array([[ 10., -31., 30.],
[ 1., 0., 0.],
[ 0., 1., 0.]])
coshm(A)
Compute the hyperbolic matrix cosine.
This routine uses expm to compute the matrix exponentials.
Parameters
----------
A : array, shape(M,M)
Returns
-------
coshA : array, shape(M,M)
Hyperbolic matrix cosine of A
cosm(A)
Compute the matrix cosine.
This routine uses expm to compute the matrix exponentials.
Parameters
----------
A : array, shape(M,M)
Returns
-------
cosA : array, shape(M,M)
Matrix cosine of A
det(a, overwrite_a=False)
Compute the determinant of a matrix
The determinant of a square matrix is a value derived arithmetically
from the coefficients of the matrix.
The determinant for a 3x3 matrix, for example, is computed as follows::
a b c
d e f = A
g h i
det(A) = a*e*i +b*f*g + c*d*h - c*e*g - b*d*i - a*f*h
Parameters
----------
a : array_like, shape (M, M)
A square matrix.
overwrite_a : bool
Allow overwriting data in a (may enhance performance).
Returns
-------
det : float or complex
Determinant of `a`.
Notes
-----
The determinant is computed via LU factorization, LAPACK routine z/dgetrf.
Examples
--------
>>> a = np.array([[1,2,3],[4,5,6],[7,8,9]])
>>> linalg.det(a)
0.0
>>> a = np.array([[0,2,3],[4,5,6],[7,8,9]])
>>> linalg.det(a)
3.0
diagsvd(s, M, N)
Construct the sigma matrix in SVD from singular values and size M, N.
Parameters
----------
s : array_like, shape (M,) or (N,)
Singular values
M : int
Size of the matrix whose singular values are `s`.
N : int
Size of the matrix whose singular values are `s`.
Returns
-------
S : array, shape (M, N)
The S-matrix in the singular value decomposition
eig(a, b=None, left=False, right=True, overwrite_a=False, overwrite_b=False)
Solve an ordinary or generalized eigenvalue problem of a square matrix.
Find eigenvalues w and right or left eigenvectors of a general matrix::
a vr[:,i] = w[i] b vr[:,i]
a.H vl[:,i] = w[i].conj() b.H vl[:,i]
where ``.H`` is the Hermitian conjugation.
Parameters
----------
a : array_like, shape (M, M)
A complex or real matrix whose eigenvalues and eigenvectors
will be computed.
b : array_like, shape (M, M), optional
Right-hand side matrix in a generalized eigenvalue problem.
Default is None, identity matrix is assumed.
left : bool, optional
Whether to calculate and return left eigenvectors. Default is False.
right : bool, optional
Whether to calculate and return right eigenvectors. Default is True.
overwrite_a : bool, optional
Whether to overwrite `a`; may improve performance. Default is False.
overwrite_b : bool, optional
Whether to overwrite `b`; may improve performance. Default is False.
Returns
-------
w : double or complex ndarray
The eigenvalues, each repeated according to its multiplicity.
Of shape (M,).
vl : double or complex ndarray
The normalized left eigenvector corresponding to the eigenvalue
``w[i]`` is the column v[:,i]. Only returned if ``left=True``.
Of shape ``(M, M)``.
vr : double or complex array
The normalized right eigenvector corresponding to the eigenvalue
``w[i]`` is the column ``vr[:,i]``. Only returned if ``right=True``.
Of shape ``(M, M)``.
Raises
------
LinAlgError
If eigenvalue computation does not converge.
See Also
--------
eigh : Eigenvalues and right eigenvectors for symmetric/Hermitian arrays.
eig_banded(a_band, lower=False, eigvals_only=False, overwrite_a_band=False, select='a', select_range=None, max_ev=0)
Solve real symmetric or complex hermitian band matrix eigenvalue problem.
Find eigenvalues w and optionally right eigenvectors v of a::
a v[:,i] = w[i] v[:,i]
v.H v = identity
The matrix a is stored in a_band either in lower diagonal or upper
diagonal ordered form:
a_band[u + i - j, j] == a[i,j] (if upper form; i <= j)
a_band[ i - j, j] == a[i,j] (if lower form; i >= j)
where u is the number of bands above the diagonal.
Example of a_band (shape of a is (6,6), u=2)::
upper form:
* * a02 a13 a24 a35
* a01 a12 a23 a34 a45
a00 a11 a22 a33 a44 a55
lower form:
a00 a11 a22 a33 a44 a55
a10 a21 a32 a43 a54 *
a20 a31 a42 a53 * *
Cells marked with * are not used.
Parameters
----------
a_band : array, shape (u+1, M)
The bands of the M by M matrix a.
lower : boolean
Is the matrix in the lower form. (Default is upper form)
eigvals_only : boolean
Compute only the eigenvalues and no eigenvectors.
(Default: calculate also eigenvectors)
overwrite_a_band:
Discard data in a_band (may enhance performance)
select: {'a', 'v', 'i'}
Which eigenvalues to calculate
====== ========================================
select calculated
====== ========================================
'a' All eigenvalues
'v' Eigenvalues in the interval (min, max]
'i' Eigenvalues with indices min <= i <= max
====== ========================================
select_range : (min, max)
Range of selected eigenvalues
max_ev : integer
For select=='v', maximum number of eigenvalues expected.
For other values of select, has no meaning.
In doubt, leave this parameter untouched.
Returns
-------
w : array, shape (M,)
The eigenvalues, in ascending order, each repeated according to its
multiplicity.
v : double or complex double array, shape (M, M)
The normalized eigenvector corresponding to the eigenvalue w[i] is
the column v[:,i].
Raises LinAlgError if eigenvalue computation does not converge
eigh(a, b=None, lower=True, eigvals_only=False, overwrite_a=False, overwrite_b=False, turbo=True, eigvals=None, type=1)
Solve an ordinary or generalized eigenvalue problem for a complex
Hermitian or real symmetric matrix.
Find eigenvalues w and optionally eigenvectors v of matrix a, where
b is positive definite::
a v[:,i] = w[i] b v[:,i]
v[i,:].conj() a v[:,i] = w[i]
v[i,:].conj() b v[:,i] = 1
Parameters
----------
a : array, shape (M, M)
A complex Hermitian or real symmetric matrix whose eigenvalues and
eigenvectors will be computed.
b : array, shape (M, M)
A complex Hermitian or real symmetric definite positive matrix in.
If omitted, identity matrix is assumed.
lower : boolean
Whether the pertinent array data is taken from the lower or upper
triangle of a. (Default: lower)
eigvals_only : boolean
Whether to calculate only eigenvalues and no eigenvectors.
(Default: both are calculated)
turbo : boolean
Use divide and conquer algorithm (faster but expensive in memory,
only for generalized eigenvalue problem and if eigvals=None)
eigvals : tuple (lo, hi)
Indexes of the smallest and largest (in ascending order) eigenvalues
and corresponding eigenvectors to be returned: 0 <= lo < hi <= M-1.
If omitted, all eigenvalues and eigenvectors are returned.
type: integer
Specifies the problem type to be solved:
type = 1: a v[:,i] = w[i] b v[:,i]
type = 2: a b v[:,i] = w[i] v[:,i]
type = 3: b a v[:,i] = w[i] v[:,i]
overwrite_a : boolean
Whether to overwrite data in a (may improve performance)
overwrite_b : boolean
Whether to overwrite data in b (may improve performance)
Returns
-------
w : real array, shape (N,)
The N (1<=N<=M) selected eigenvalues, in ascending order, each
repeated according to its multiplicity.
(if eigvals_only == False)
v : complex array, shape (M, N)
The normalized selected eigenvector corresponding to the
eigenvalue w[i] is the column v[:,i]. Normalization:
type 1 and 3: v.conj() a v = w
type 2: inv(v).conj() a inv(v) = w
type = 1 or 2: v.conj() b v = I
type = 3 : v.conj() inv(b) v = I
Raises LinAlgError if eigenvalue computation does not converge,
an error occurred, or b matrix is not definite positive. Note that
if input matrices are not symmetric or hermitian, no error is reported
but results will be wrong.
See Also
--------
eig : eigenvalues and right eigenvectors for non-symmetric arrays
eigvals(a, b=None, overwrite_a=False)
Compute eigenvalues from an ordinary or generalized eigenvalue problem.
Find eigenvalues of a general matrix::
a vr[:,i] = w[i] b vr[:,i]
Parameters
----------
a : array_like, shape (M, M)
A complex or real matrix whose eigenvalues and eigenvectors
will be computed.
b : array_like, shape (M, M), optional
Right-hand side matrix in a generalized eigenvalue problem.
If omitted, identity matrix is assumed.
overwrite_a : boolean, optional
Whether to overwrite data in a (may improve performance)
Returns
-------
w : double or complex ndarray, shape (M,)
The eigenvalues, each repeated according to its multiplicity,
but not in any specific order. Of shape (M,).
Raises
------
LinAlgError
If eigenvalue computation does not converge
See Also
--------
eigvalsh : eigenvalues of symmetric or Hermitian arrays,
eig : eigenvalues and right eigenvectors of general arrays.
eigh : eigenvalues and eigenvectors of symmetric/Hermitian arrays.
eigvals_banded(a_band, lower=False, overwrite_a_band=False, select='a', select_range=None)
Solve real symmetric or complex hermitian band matrix eigenvalue problem.
Find eigenvalues w of a::
a v[:,i] = w[i] v[:,i]
v.H v = identity
The matrix a is stored in a_band either in lower diagonal or upper
diagonal ordered form:
a_band[u + i - j, j] == a[i,j] (if upper form; i <= j)
a_band[ i - j, j] == a[i,j] (if lower form; i >= j)
where u is the number of bands above the diagonal.
Example of a_band (shape of a is (6,6), u=2)::
upper form:
* * a02 a13 a24 a35
* a01 a12 a23 a34 a45
a00 a11 a22 a33 a44 a55
lower form:
a00 a11 a22 a33 a44 a55
a10 a21 a32 a43 a54 *
a20 a31 a42 a53 * *
Cells marked with * are not used.
Parameters
----------
a_band : array, shape (u+1, M)
The bands of the M by M matrix a.
lower : boolean
Is the matrix in the lower form. (Default is upper form)
overwrite_a_band:
Discard data in a_band (may enhance performance)
select: {'a', 'v', 'i'}
Which eigenvalues to calculate
====== ========================================
select calculated
====== ========================================
'a' All eigenvalues
'v' Eigenvalues in the interval (min, max]
'i' Eigenvalues with indices min <= i <= max
====== ========================================
select_range : (min, max)
Range of selected eigenvalues
Returns
-------
w : array, shape (M,)
The eigenvalues, in ascending order, each repeated according to its
multiplicity.
Raises LinAlgError if eigenvalue computation does not converge
See Also
--------
eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian band matrices
eigvals : eigenvalues of general arrays
eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
eig : eigenvalues and right eigenvectors for non-symmetric arrays
eigvalsh(a, b=None, lower=True, overwrite_a=False, overwrite_b=False, turbo=True, eigvals=None, type=1)
Solve an ordinary or generalized eigenvalue problem for a complex
Hermitian or real symmetric matrix.
Find eigenvalues w of matrix a, where b is positive definite::
a v[:,i] = w[i] b v[:,i]
v[i,:].conj() a v[:,i] = w[i]
v[i,:].conj() b v[:,i] = 1
Parameters
----------
a : array, shape (M, M)
A complex Hermitian or real symmetric matrix whose eigenvalues and
eigenvectors will be computed.
b : array, shape (M, M)
A complex Hermitian or real symmetric definite positive matrix in.
If omitted, identity matrix is assumed.
lower : boolean
Whether the pertinent array data is taken from the lower or upper
triangle of a. (Default: lower)
turbo : boolean
Use divide and conquer algorithm (faster but expensive in memory,
only for generalized eigenvalue problem and if eigvals=None)
eigvals : tuple (lo, hi)
Indexes of the smallest and largest (in ascending order) eigenvalues
and corresponding eigenvectors to be returned: 0 <= lo < hi <= M-1.
If omitted, all eigenvalues and eigenvectors are returned.
type: integer
Specifies the problem type to be solved:
type = 1: a v[:,i] = w[i] b v[:,i]
type = 2: a b v[:,i] = w[i] v[:,i]
type = 3: b a v[:,i] = w[i] v[:,i]
overwrite_a : boolean
Whether to overwrite data in a (may improve performance)
overwrite_b : boolean
Whether to overwrite data in b (may improve performance)
Returns
-------
w : real array, shape (N,)
The N (1<=N<=M) selected eigenvalues, in ascending order, each
repeated according to its multiplicity.
Raises LinAlgError if eigenvalue computation does not converge,
an error occurred, or b matrix is not definite positive. Note that
if input matrices are not symmetric or hermitian, no error is reported
but results will be wrong.
See Also
--------
eigvals : eigenvalues of general arrays
eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
eig : eigenvalues and right eigenvectors for non-symmetric arrays
expm(A, q=False)
Compute the matrix exponential using Pade approximation.
Parameters
----------
A : array, shape(M,M)
Matrix to be exponentiated
Returns
-------
expA : array, shape(M,M)
Matrix exponential of A
References
----------
N. J. Higham,
"The Scaling and Squaring Method for the Matrix Exponential Revisited",
SIAM. J. Matrix Anal. & Appl. 26, 1179 (2005).
expm2(A)
Compute the matrix exponential using eigenvalue decomposition.
Parameters
----------
A : array, shape(M,M)
Matrix to be exponentiated
Returns
-------
expA : array, shape(M,M)
Matrix exponential of A
expm3(A, q=20)
Compute the matrix exponential using Taylor series.
Parameters
----------
A : array, shape(M,M)
Matrix to be exponentiated
q : integer
Order of the Taylor series
Returns
-------
expA : array, shape(M,M)
Matrix exponential of A
funm(A, func, disp=True)
Evaluate a matrix function specified by a callable.
Returns the value of matrix-valued function f at A. The function f
is an extension of the scalar-valued function func to matrices.
Parameters
----------
A : array, shape(M,M)
Matrix at which to evaluate the function
func : callable
Callable object that evaluates a scalar function f.
Must be vectorized (eg. using vectorize).
disp : boolean
Print warning if error in the result is estimated large
instead of returning estimated error. (Default: True)
Returns
-------
fA : array, shape(M,M)
Value of the matrix function specified by func evaluated at A
(if disp == False)
errest : float
1-norm of the estimated error, ||err||_1 / ||A||_1
get_blas_funcs(names, arrays=(), dtype=None)
Return available BLAS function objects from names.
Arrays are used to determine the optimal prefix of BLAS routines.
Parameters
----------
names : str or sequence of str
Name(s) of BLAS functions withouth type prefix.
arrays : sequency of ndarrays, optional
Arrays can be given to determine optiomal prefix of BLAS
routines. If not given, double-precision routines will be
used, otherwise the most generic type in arrays will be used.
dtype : str or dtype, optional
Data-type specifier. Not used if `arrays` is non-empty.
Returns
-------
funcs : list
List containing the found function(s).
Notes
-----
This routines automatically chooses between Fortran/C
interfaces. Fortran code is used whenever possible for arrays with
column major order. In all other cases, C code is preferred.
In BLAS, the naming convention is that all functions start with a
type prefix, which depends on the type of the principal
matrix. These can be one of {'s', 'd', 'c', 'z'} for the numpy
types {float32, float64, complex64, complex128} respectevely, and
are stored in attribute `typecode` of the returned functions.
hadamard(n, dtype=<type 'int'>)
Construct a Hadamard matrix.
`hadamard(n)` constructs an n-by-n Hadamard matrix, using Sylvester's
construction. `n` must be a power of 2.
Parameters
----------
n : int
The order of the matrix. `n` must be a power of 2.
dtype : numpy dtype
The data type of the array to be constructed.
Returns
-------
H : ndarray with shape (n, n)
The Hadamard matrix.
Notes
-----
.. versionadded:: 0.8.0
Examples
--------
>>> from scipy.linalg import hadamard
>>> hadamard(2, dtype=complex)
array([[ 1.+0.j, 1.+0.j],
[ 1.+0.j, -1.-0.j]])
>>> hadamard(4)
array([[ 1, 1, 1, 1],
[ 1, -1, 1, -1],
[ 1, 1, -1, -1],
[ 1, -1, -1, 1]])
hankel(c, r=None)
Construct a Hankel matrix.
The Hankel matrix has constant anti-diagonals, with `c` as its
first column and `r` as its last row. If `r` is not given, then
`r = zeros_like(c)` is assumed.
Parameters
----------
c : array_like
First column of the matrix. Whatever the actual shape of `c`, it
will be converted to a 1-D array.
r : array_like, 1D
Last row of the matrix. If None, ``r = zeros_like(c)`` is assumed.
r[0] is ignored; the last row of the returned matrix is
``[c[-1], r[1:]]``. Whatever the actual shape of `r`, it will be
converted to a 1-D array.
Returns
-------
A : array, shape (len(c), len(r))
The Hankel matrix. Dtype is the same as ``(c[0] + r[0]).dtype``.
See also
--------
toeplitz : Toeplitz matrix
circulant : circulant matrix
Examples
--------
>>> from scipy.linalg import hankel
>>> hankel([1, 17, 99])
array([[ 1, 17, 99],
[17, 99, 0],
[99, 0, 0]])
>>> hankel([1,2,3,4], [4,7,7,8,9])
array([[1, 2, 3, 4, 7],
[2, 3, 4, 7, 7],
[3, 4, 7, 7, 8],
[4, 7, 7, 8, 9]])
hessenberg(a, calc_q=False, overwrite_a=False)
Compute Hessenberg form of a matrix.
The Hessenberg decomposition is::
A = Q H Q^H
where `Q` is unitary/orthogonal and `H` has only zero elements below
the first sub-diagonal.
Parameters
----------
a : ndarray
Matrix to bring into Hessenberg form, of shape ``(M,M)``.
calc_q : bool, optional
Whether to compute the transformation matrix. Default is False.
overwrite_a : bool, optional
Whether to overwrite `a`; may improve performance.
Default is False.
Returns
-------
H : ndarray
Hessenberg form of `a`, of shape (M,M).
Q : ndarray
Unitary/orthogonal similarity transformation matrix ``A = Q H Q^H``.
Only returned if ``calc_q=True``. Of shape (M,M).
hilbert(n)
Create a Hilbert matrix of order `n`.
Returns the `n` by `n` array with entries `h[i,j] = 1 / (i + j + 1)`.
Parameters
----------
n : int
The size of the array to create.
Returns
-------
h : ndarray with shape (n, n)
The Hilbert matrix.
See Also
--------
invhilbert : Compute the inverse of a Hilbert matrix.
Notes
-----
.. versionadded:: 0.10.0
Examples
--------
>>> from scipy.linalg import hilbert
>>> hilbert(3)
array([[ 1. , 0.5 , 0.33333333],
[ 0.5 , 0.33333333, 0.25 ],
[ 0.33333333, 0.25 , 0.2 ]])
inv(a, overwrite_a=False)
Compute the inverse of a matrix.
Parameters
----------
a : array_like
Square matrix to be inverted.
overwrite_a : bool, optional
Discard data in `a` (may improve performance). Default is False.
Returns
-------
ainv : ndarray
Inverse of the matrix `a`.
Raises
------
LinAlgError :
If `a` is singular.
ValueError :
If `a` is not square, or not 2-dimensional.
Examples
--------
>>> a = np.array([[1., 2.], [3., 4.]])
>>> sp.linalg.inv(a)
array([[-2. , 1. ],
[ 1.5, -0.5]])
>>> np.dot(a, sp.linalg.inv(a))
array([[ 1., 0.],
[ 0., 1.]])
invhilbert(n, exact=False)
Compute the inverse of the Hilbert matrix of order `n`.
The entries in the inverse of a Hilbert matrix are integers. When `n`
is greater than 14, some entries in the inverse exceed the upper limit
of 64 bit integers. The `exact` argument provides two options for
dealing with these large integers.
Parameters
----------
n : int
The order of the Hilbert matrix.
exact : bool
If False, the data type of the array that is returned is np.float64,
and the array is an approximation of the inverse.
If True, the array is the exact integer inverse array. To represent
the exact inverse when n > 14, the returned array is an object array
of long integers. For n <= 14, the exact inverse is returned as an
array with data type np.int64.
Returns
-------
invh : ndarray with shape (n, n)
The data type of the array is np.float64 if `exact` is False.
If `exact` is True, the data type is either np.int64 (for n <= 14)
or object (for n > 14). In the latter case, the objects in the
array will be long integers.
See Also
--------
hilbert : Create a Hilbert matrix.
Notes
-----
.. versionadded:: 0.10.0
Examples
--------
>>> from scipy.linalg import invhilbert
>>> invhilbert(4)
array([[ 16., -120., 240., -140.],
[ -120., 1200., -2700., 1680.],
[ 240., -2700., 6480., -4200.],
[ -140., 1680., -4200., 2800.]])
>>> invhilbert(4, exact=True)
array([[ 16, -120, 240, -140],
[ -120, 1200, -2700, 1680],
[ 240, -2700, 6480, -4200],
[ -140, 1680, -4200, 2800]], dtype=int64)
>>> invhilbert(16)[7,7]
4.2475099528537506e+19
>>> invhilbert(16, exact=True)[7,7]
42475099528537378560L
kron(a, b)
Kronecker product of a and b.
The result is the block matrix::
a[0,0]*b a[0,1]*b ... a[0,-1]*b
a[1,0]*b a[1,1]*b ... a[1,-1]*b
...
a[-1,0]*b a[-1,1]*b ... a[-1,-1]*b
Parameters
----------
a : array, shape (M, N)
b : array, shape (P, Q)
Returns
-------
A : array, shape (M*P, N*Q)
Kronecker product of a and b
Examples
--------
>>> from numpy import array
>>> from scipy.linalg import kron
>>> kron(array([[1,2],[3,4]]), array([[1,1,1]]))
array([[1, 1, 1, 2, 2, 2],
[3, 3, 3, 4, 4, 4]])
leslie(f, s)
Create a Leslie matrix.
Given the length n array of fecundity coefficients `f` and the length
n-1 array of survival coefficents `s`, return the associated Leslie matrix.
Parameters
----------
f : array_like
The "fecundity" coefficients, has to be 1-D.
s : array_like
The "survival" coefficients, has to be 1-D. The length of `s`
must be one less than the length of `f`, and it must be at least 1.
Returns
-------
L : ndarray
Returns a 2-D ndarray of shape ``(n, n)``, where `n` is the
length of `f`. The array is zero except for the first row,
which is `f`, and the first sub-diagonal, which is `s`.
The data-type of the array will be the data-type of ``f[0]+s[0]``.
Notes
-----
.. versionadded:: 0.8.0
The Leslie matrix is used to model discrete-time, age-structured
population growth [1]_ [2]_. In a population with `n` age classes, two sets
of parameters define a Leslie matrix: the `n` "fecundity coefficients",
which give the number of offspring per-capita produced by each age
class, and the `n` - 1 "survival coefficients", which give the
per-capita survival rate of each age class.
References
----------
.. [1] P. H. Leslie, On the use of matrices in certain population
mathematics, Biometrika, Vol. 33, No. 3, 183--212 (Nov. 1945)
.. [2] P. H. Leslie, Some further notes on the use of matrices in
population mathematics, Biometrika, Vol. 35, No. 3/4, 213--245
(Dec. 1948)
Examples
--------
>>> from scipy.linalg import leslie
>>> leslie([0.1, 2.0, 1.0, 0.1], [0.2, 0.8, 0.7])
array([[ 0.1, 2. , 1. , 0.1],
[ 0.2, 0. , 0. , 0. ],
[ 0. , 0.8, 0. , 0. ],
[ 0. , 0. , 0.7, 0. ]])
logm(A, disp=True)
Compute matrix logarithm.
The matrix logarithm is the inverse of expm: expm(logm(A)) == A
Parameters
----------
A : array, shape(M,M)
Matrix whose logarithm to evaluate
disp : boolean
Print warning if error in the result is estimated large
instead of returning estimated error. (Default: True)
Returns
-------
logA : array, shape(M,M)
Matrix logarithm of A
(if disp == False)
errest : float
1-norm of the estimated error, ||err||_1 / ||A||_1
lstsq(a, b, cond=None, overwrite_a=False, overwrite_b=False)
Compute least-squares solution to equation Ax = b.
Compute a vector x such that the 2-norm ``|b - A x|`` is minimized.
Parameters
----------
a : array, shape (M, N)
Left hand side matrix (2-D array).
b : array, shape (M,) or (M, K)
Right hand side matrix or vector (1-D or 2-D array).
cond : float, optional
Cutoff for 'small' singular values; used to determine effective
rank of a. Singular values smaller than
``rcond * largest_singular_value`` are considered zero.
overwrite_a : bool, optional
Discard data in `a` (may enhance performance). Default is False.
overwrite_b : bool, optional
Discard data in `b` (may enhance performance). Default is False.
Returns
-------
x : array, shape (N,) or (N, K) depending on shape of b
Least-squares solution.
residues : ndarray, shape () or (1,) or (K,)
Sums of residues, squared 2-norm for each column in ``b - a x``.
If rank of matrix a is < N or > M this is an empty array.
If b was 1-D, this is an (1,) shape array, otherwise the shape is (K,).
rank : int
Effective rank of matrix `a`.
s : array, shape (min(M,N),)
Singular values of `a`. The condition number of a is
``abs(s[0]/s[-1])``.
Raises
------
LinAlgError :
If computation does not converge.
See Also
--------
optimize.nnls : linear least squares with non-negativity constraint
lu(a, permute_l=False, overwrite_a=False)
Compute pivoted LU decompostion of a matrix.
The decomposition is::
A = P L U
where P is a permutation matrix, L lower triangular with unit
diagonal elements, and U upper triangular.
Parameters
----------
a : array, shape (M, N)
Array to decompose
permute_l : boolean
Perform the multiplication P*L (Default: do not permute)
overwrite_a : boolean
Whether to overwrite data in a (may improve performance)
Returns
-------
(If permute_l == False)
p : array, shape (M, M)
Permutation matrix
l : array, shape (M, K)
Lower triangular or trapezoidal matrix with unit diagonal.
K = min(M, N)
u : array, shape (K, N)
Upper triangular or trapezoidal matrix
(If permute_l == True)
pl : array, shape (M, K)
Permuted L matrix.
K = min(M, N)
u : array, shape (K, N)
Upper triangular or trapezoidal matrix
Notes
-----
This is a LU factorization routine written for Scipy.
lu_factor(a, overwrite_a=False)
Compute pivoted LU decomposition of a matrix.
The decomposition is::
A = P L U
where P is a permutation matrix, L lower triangular with unit
diagonal elements, and U upper triangular.
Parameters
----------
a : array, shape (M, M)
Matrix to decompose
overwrite_a : boolean
Whether to overwrite data in A (may increase performance)
Returns
-------
lu : array, shape (N, N)
Matrix containing U in its upper triangle, and L in its lower triangle.
The unit diagonal elements of L are not stored.
piv : array, shape (N,)
Pivot indices representing the permutation matrix P:
row i of matrix was interchanged with row piv[i].
See also
--------
lu_solve : solve an equation system using the LU factorization of a matrix
Notes
-----
This is a wrapper to the ``*GETRF`` routines from LAPACK.
lu_solve((lu, piv), b, trans=0, overwrite_b=False)
Solve an equation system, a x = b, given the LU factorization of a
Parameters
----------
(lu, piv)
Factorization of the coefficient matrix a, as given by lu_factor
b : array
Right-hand side
trans : {0, 1, 2}
Type of system to solve:
===== =========
trans system
===== =========
0 a x = b
1 a^T x = b
2 a^H x = b
===== =========
Returns
-------
x : array
Solution to the system
See also
--------
lu_factor : LU factorize a matrix
norm(a, ord=None)
Matrix or vector norm.
This function is able to return one of seven different matrix norms,
or one of an infinite number of vector norms (described below), depending
on the value of the ``ord`` parameter.
Parameters
----------
x : {(M,), (M, N)} array_like
Input array.
ord : {non-zero int, inf, -inf, 'fro'}, optional
Order of the norm (see table under ``Notes``). inf means numpy's
`inf` object.
Returns
-------
n : float
Norm of the matrix or vector.
Notes
-----
For values of ``ord <= 0``, the result is, strictly speaking, not a
mathematical 'norm', but it may still be useful for various numerical
purposes.
The following norms can be calculated:
===== ============================ ==========================
ord norm for matrices norm for vectors
===== ============================ ==========================
None Frobenius norm 2-norm
'fro' Frobenius norm --
inf max(sum(abs(x), axis=1)) max(abs(x))
-inf min(sum(abs(x), axis=1)) min(abs(x))
0 -- sum(x != 0)
1 max(sum(abs(x), axis=0)) as below
-1 min(sum(abs(x), axis=0)) as below
2 2-norm (largest sing. value) as below
-2 smallest singular value as below
other -- sum(abs(x)**ord)**(1./ord)
===== ============================ ==========================
The Frobenius norm is given by [1]_:
:math:`||A||_F = [\sum_{i,j} abs(a_{i,j})^2]^{1/2}`
References
----------
.. [1] G. H. Golub and C. F. Van Loan, *Matrix Computations*,
Baltimore, MD, Johns Hopkins University Press, 1985, pg. 15
Examples
--------
>>> from numpy import linalg as LA
>>> a = np.arange(9) - 4
>>> a
array([-4, -3, -2, -1, 0, 1, 2, 3, 4])
>>> b = a.reshape((3, 3))
>>> b
array([[-4, -3, -2],
[-1, 0, 1],
[ 2, 3, 4]])
>>> LA.norm(a)
7.745966692414834
>>> LA.norm(b)
7.745966692414834
>>> LA.norm(b, 'fro')
7.745966692414834
>>> LA.norm(a, np.inf)
4
>>> LA.norm(b, np.inf)
9
>>> LA.norm(a, -np.inf)
0
>>> LA.norm(b, -np.inf)
2
>>> LA.norm(a, 1)
20
>>> LA.norm(b, 1)
7
>>> LA.norm(a, -1)
-4.6566128774142013e-010
>>> LA.norm(b, -1)
6
>>> LA.norm(a, 2)
7.745966692414834
>>> LA.norm(b, 2)
7.3484692283495345
>>> LA.norm(a, -2)
nan
>>> LA.norm(b, -2)
1.8570331885190563e-016
>>> LA.norm(a, 3)
5.8480354764257312
>>> LA.norm(a, -3)
nan
orth(A)
Construct an orthonormal basis for the range of A using SVD
Parameters
----------
A : array, shape (M, N)
Returns
-------
Q : array, shape (M, K)
Orthonormal basis for the range of A.
K = effective rank of A, as determined by automatic cutoff
See also
--------
svd : Singular value decomposition of a matrix
pascal(n, kind='symmetric', exact=True)
Returns the n x n Pascal matrix.
The Pascal matrix is a matrix containing the binomial coefficients as
its elements.
Parameters
----------
n : int
The size of the matrix to create; that is, the result is an n x n
matrix.
kind : str, optional
Must be one of 'symmetric', 'lower', or 'upper'.
Default is 'symmetric'.
exact : bool, optional
If `exact` is True, the result is either an array of type
numpy.uint64 (if n <= 35) or an object array of Python long integers.
If `exact` is False, the coefficients in the matrix are computed using
`scipy.misc.comb` with `exact=False`. The result will be a floating
point array, and the values in the array will not be the exact
coefficients, but this version is much faster than `exact=True`.
Returns
-------
p : 2-D ndarray
The Pascal matrix.
Notes
-----
.. versionadded:: 0.11.0
See http://en.wikipedia.org/wiki/Pascal_matrix for more information
about Pascal matrices.
Examples
--------
>>> from scipy.linalg import pascal
>>> pascal(4)
array([[ 1, 1, 1, 1],
[ 1, 2, 3, 4],
[ 1, 3, 6, 10],
[ 1, 4, 10, 20]], dtype=uint64)
>>> pascal(4, kind='lower')
array([[1, 0, 0, 0],
[1, 1, 0, 0],
[1, 2, 1, 0],
[1, 3, 3, 1]], dtype=uint64)
>>> pascal(50)[-1, -1]
25477612258980856902730428600L
>>> from scipy.misc import comb
>>> comb(98, 49, exact=True)
25477612258980856902730428600L
pinv(a, cond=None, rcond=None)
Compute the (Moore-Penrose) pseudo-inverse of a matrix.
Calculate a generalized inverse of a matrix using a least-squares
solver.
Parameters
----------
a : array, shape (M, N)
Matrix to be pseudo-inverted.
cond, rcond : float, optional
Cutoff for 'small' singular values in the least-squares solver.
Singular values smaller than ``rcond * largest_singular_value``
are considered zero.
Returns
-------
B : array, shape (N, M)
The pseudo-inverse of matrix `a`.
Raises
------
LinAlgError
If computation does not converge.
Examples
--------
>>> a = np.random.randn(9, 6)
>>> B = linalg.pinv(a)
>>> np.allclose(a, dot(a, dot(B, a)))
True
>>> np.allclose(B, dot(B, dot(a, B)))
True
pinv2(a, cond=None, rcond=None)
Compute the (Moore-Penrose) pseudo-inverse of a matrix.
Calculate a generalized inverse of a matrix using its
singular-value decomposition and including all 'large' singular
values.
Parameters
----------
a : array, shape (M, N)
Matrix to be pseudo-inverted.
cond, rcond : float or None
Cutoff for 'small' singular values.
Singular values smaller than ``rcond*largest_singular_value``
are considered zero.
If None or -1, suitable machine precision is used.
Returns
-------
B : array, shape (N, M)
The pseudo-inverse of matrix `a`.
Raises
------
LinAlgError
If SVD computation does not converge.
Examples
--------
>>> a = np.random.randn(9, 6)
>>> B = linalg.pinv2(a)
>>> np.allclose(a, dot(a, dot(B, a)))
True
>>> np.allclose(B, dot(B, dot(a, B)))
True
qr(a, overwrite_a=False, lwork=None, mode='full', pivoting=False)
Compute QR decomposition of a matrix.
Calculate the decomposition ``A = Q R`` where Q is unitary/orthogonal
and R upper triangular.
Parameters
----------
a : array, shape (M, N)
Matrix to be decomposed
overwrite_a : bool, optional
Whether data in a is overwritten (may improve performance)
lwork : int, optional
Work array size, lwork >= a.shape[1]. If None or -1, an optimal size
is computed.
mode : {'full', 'r', 'economic', 'raw'}
Determines what information is to be returned: either both Q and R
('full', default), only R ('r') or both Q and R but computed in
economy-size ('economic', see Notes). The final option 'raw'
(added in Scipy 0.11) makes the function return two matrixes
(Q, TAU) in the internal format used by LAPACK.
pivoting : bool, optional
Whether or not factorization should include pivoting for rank-revealing
qr decomposition. If pivoting, compute the decomposition
``A P = Q R`` as above, but where P is chosen such that the diagonal
of R is non-increasing.
Returns
-------
Q : float or complex ndarray
Of shape (M, M), or (M, K) for ``mode='economic'``. Not returned
if ``mode='r'``.
R : float or complex ndarray
Of shape (M, N), or (K, N) for ``mode='economic'``. ``K = min(M, N)``.
P : integer ndarray
Of shape (N,) for ``pivoting=True``. Not returned if
``pivoting=False``.
Raises
------
LinAlgError
Raised if decomposition fails
Notes
-----
This is an interface to the LAPACK routines dgeqrf, zgeqrf,
dorgqr, zungqr, dgeqp3, and zgeqp3.
If ``mode=economic``, the shapes of Q and R are (M, K) and (K, N) instead
of (M,M) and (M,N), with ``K=min(M,N)``.
Examples
--------
>>> from scipy import random, linalg, dot, diag, all, allclose
>>> a = random.randn(9, 6)
>>> q, r = linalg.qr(a)
>>> allclose(a, np.dot(q, r))
True
>>> q.shape, r.shape
((9, 9), (9, 6))
>>> r2 = linalg.qr(a, mode='r')
>>> allclose(r, r2)
True
>>> q3, r3 = linalg.qr(a, mode='economic')
>>> q3.shape, r3.shape
((9, 6), (6, 6))
>>> q4, r4, p4 = linalg.qr(a, pivoting=True)
>>> d = abs(diag(r4))
>>> all(d[1:] <= d[:-1])
True
>>> allclose(a[:, p4], dot(q4, r4))
True
>>> q4.shape, r4.shape, p4.shape
((9, 9), (9, 6), (6,))
>>> q5, r5, p5 = linalg.qr(a, mode='economic', pivoting=True)
>>> q5.shape, r5.shape, p5.shape
((9, 6), (6, 6), (6,))
qr_multiply(a, c, mode='right', pivoting=False, conjugate=False, overwrite_a=False, overwrite_c=False)
Calculate the QR decomposition and multiply Q with a matrix.
Calculate the decomposition ``A = Q R`` where Q is unitary/orthogonal
and R upper triangular. Multiply Q with a vector or a matrix c.
.. versionadded:: 0.11
Parameters
----------
a : ndarray, shape (M, N)
Matrix to be decomposed
c : ndarray, one- or two-dimensional
calculate the product of c and q, depending on the mode:
mode : {'left', 'right'}
``dot(Q, c)`` is returned if mode is 'left',
``dot(c, Q)`` is returned if mode is 'right'.
The shape of c must be appropriate for the matrix multiplications,
if mode is 'left', ``min(a.shape) == c.shape[0]``,
if mode is 'right', ``a.shape[0] == c.shape[1]``.
pivoting : bool, optional
Whether or not factorization should include pivoting for rank-revealing
qr decomposition, see the documentation of qr.
conjugate : bool, optional
Whether Q should be complex-conjugated. This might be faster
than explicit conjugation.
overwrite_a : bool, optional
Whether data in a is overwritten (may improve performance)
overwrite_c: bool, optional
Whether data in c is overwritten (may improve performance).
If this is used, c must be big enough to keep the result,
i.e. c.shape[0] = a.shape[0] if mode is 'left'.
Returns
-------
CQ : float or complex ndarray
the product of Q and c, as defined in mode
R : float or complex ndarray
Of shape (K, N), ``K = min(M, N)``.
P : ndarray of ints
Of shape (N,) for ``pivoting=True``. Not returned if ``pivoting=False``.
Raises
------
LinAlgError
Raised if decomposition fails
Notes
-----
This is an interface to the LAPACK routines dgeqrf, zgeqrf,
dormqr, zunmqr, dgeqp3, and zgeqp3.
qr_old(*args, **kwds)
`qr_old` is deprecated!
Compute QR decomposition of a matrix.
Calculate the decomposition :lm:`A = Q R` where Q is unitary/orthogonal
and R upper triangular.
Parameters
----------
a : array, shape (M, N)
Matrix to be decomposed
overwrite_a : boolean
Whether data in a is overwritten (may improve performance)
lwork : integer
Work array size, lwork >= a.shape[1]. If None or -1, an optimal size
is computed.
Returns
-------
Q : float or complex array, shape (M, M)
R : float or complex array, shape (M, N)
Size K = min(M, N)
Raises LinAlgError if decomposition fails
qz(A, B, output='real', lwork=None, sort=None, overwrite_a=False, overwrite_b=False)
QZ decompostion for generalized eigenvalues of a pair of matrices.
The QZ, or generalized Schur, decomposition for a pair of N x N
nonsymmetric matrices (A,B) is::
(A,B) = (Q*AA*Z', Q*BB*Z')
where AA, BB is in generalized Schur form if BB is upper-triangular
with non-negative diagonal and AA is upper-triangular, or for real QZ
decomposition (``output='real'``) block upper triangular with 1x1
and 2x2 blocks. In this case, the 1x1 blocks correspond to real
generalized eigenvalues and 2x2 blocks are 'standardized' by making
the corresponding elements of BB have the form::
[ a 0 ]
[ 0 b ]
and the pair of corresponding 2x2 blocks in AA and BB will have a complex
conjugate pair of generalized eigenvalues. If (``output='complex'``) or
A and B are complex matrices, Z' denotes the conjugate-transpose of Z.
Q and Z are unitary matrices.
Parameters
----------
A : array_like, shape (N,N)
2-D array to decompose.
B : array_like, shape (N,N)
2-D array to decompose.
output : {'real','complex'}, optional
Construct the real or complex QZ decomposition for real matrices.
Default is 'real'.
lwork : int, optional
Work array size. If None or -1, it is automatically computed.
sort : {None, callable, 'lhp', 'rhp', 'iuc', 'ouc'}, optional
NOTE: THIS INPUT IS DISABLED FOR NOW, IT DOESN'T WORK WELL ON WINDOWS.
Specifies whether the upper eigenvalues should be sorted. A callable
may be passed that, given a eigenvalue, returns a boolean denoting
whether the eigenvalue should be sorted to the top-left (True). For
real matrix pairs, the sort function takes three real arguments
(alphar, alphai, beta). The eigenvalue x = (alphar + alphai*1j)/beta.
For complex matrix pairs or output='complex', the sort function
takes two complex arguments (alpha, beta). The eigenvalue
x = (alpha/beta).
Alternatively, string parameters may be used:
- 'lhp' Left-hand plane (x.real < 0.0)
- 'rhp' Right-hand plane (x.real > 0.0)
- 'iuc' Inside the unit circle (x*x.conjugate() <= 1.0)
- 'ouc' Outside the unit circle (x*x.conjugate() > 1.0)
Defaults to None (no sorting).
Returns
-------
AA : ndarray, shape (N,N)
Generalized Schur form of A.
BB : ndarray, shape (N,N)
Generalized Schur form of B.
Q : ndarray, shape (N,N)
The left Schur vectors.
Z : ndarray, shape (N,N)
The right Schur vectors.
sdim : int, optional
If sorting was requested, a fifth return value will contain the
number of eigenvalues for which the sort condition was True.
Notes
-----
Q is transposed versus the equivalent function in Matlab.
.. versionadded:: 0.11.0
Examples
--------
>>> from scipy import linalg
>>> np.random.seed(1234)
>>> A = np.arange(9).reshape((3, 3))
>>> B = np.random.randn(3, 3)
>>> AA, BB, Q, Z = linalg.qz(A, B)
>>> AA
array([[-13.40928183, -4.62471562, 1.09215523],
[ 0. , 0. , 1.22805978],
[ 0. , 0. , 0.31973817]])
>>> BB
array([[ 0.33362547, -1.37393632, 0.02179805],
[ 0. , 1.68144922, 0.74683866],
[ 0. , 0. , 0.9258294 ]])
>>> Q
array([[ 0.14134727, -0.97562773, 0.16784365],
[ 0.49835904, -0.07636948, -0.86360059],
[ 0.85537081, 0.20571399, 0.47541828]])
>>> Z
array([[-0.24900855, -0.51772687, 0.81850696],
[-0.79813178, 0.58842606, 0.12938478],
[-0.54861681, -0.6210585 , -0.55973739]])
rq(a, overwrite_a=False, lwork=None, mode='full')
Compute RQ decomposition of a square real matrix.
Calculate the decomposition :lm:`A = R Q` where Q is unitary/orthogonal
and R upper triangular.
Parameters
----------
a : array, shape (M, M)
Matrix to be decomposed
overwrite_a : boolean
Whether data in a is overwritten (may improve performance)
lwork : integer
Work array size, lwork >= a.shape[1]. If None or -1, an optimal size
is computed.
mode : {'full', 'r', 'economic'}
Determines what information is to be returned: either both Q and R
('full', default), only R ('r') or both Q and R but computed in
economy-size ('economic', see Notes).
Returns
-------
R : float array, shape (M, N)
Q : float or complex array, shape (M, M)
Raises LinAlgError if decomposition fails
Examples
--------
>>> from scipy import linalg
>>> from numpy import random, dot, allclose
>>> a = random.randn(6, 9)
>>> r, q = linalg.rq(a)
>>> allclose(a, dot(r, q))
True
>>> r.shape, q.shape
((6, 9), (9, 9))
>>> r2 = linalg.rq(a, mode='r')
>>> allclose(r, r2)
True
>>> r3, q3 = linalg.rq(a, mode='economic')
>>> r3.shape, q3.shape
((6, 6), (6, 9))
rsf2csf(T, Z)
Convert real Schur form to complex Schur form.
Convert a quasi-diagonal real-valued Schur form to the upper triangular
complex-valued Schur form.
Parameters
----------
T : array, shape (M, M)
Real Schur form of the original matrix
Z : array, shape (M, M)
Schur transformation matrix
Returns
-------
T : array, shape (M, M)
Complex Schur form of the original matrix
Z : array, shape (M, M)
Schur transformation matrix corresponding to the complex form
See also
--------
schur : Schur decompose a matrix
schur(a, output='real', lwork=None, overwrite_a=False, sort=None)
Compute Schur decomposition of a matrix.
The Schur decomposition is::
A = Z T Z^H
where Z is unitary and T is either upper-triangular, or for real
Schur decomposition (output='real'), quasi-upper triangular. In
the quasi-triangular form, 2x2 blocks describing complex-valued
eigenvalue pairs may extrude from the diagonal.
Parameters
----------
a : ndarray, shape (M, M)
Matrix to decompose
output : {'real', 'complex'}, optional
Construct the real or complex Schur decomposition (for real matrices).
lwork : int, optional
Work array size. If None or -1, it is automatically computed.
overwrite_a : bool, optional
Whether to overwrite data in a (may improve performance).
sort : {None, callable, 'lhp', 'rhp', 'iuc', 'ouc'}, optional
Specifies whether the upper eigenvalues should be sorted. A callable
may be passed that, given a eigenvalue, returns a boolean denoting
whether the eigenvalue should be sorted to the top-left (True).
Alternatively, string parameters may be used::
'lhp' Left-hand plane (x.real < 0.0)
'rhp' Right-hand plane (x.real > 0.0)
'iuc' Inside the unit circle (x*x.conjugate() <= 1.0)
'ouc' Outside the unit circle (x*x.conjugate() > 1.0)
Defaults to None (no sorting).
Returns
-------
T : ndarray, shape (M, M)
Schur form of A. It is real-valued for the real Schur decomposition.
Z : ndarray, shape (M, M)
An unitary Schur transformation matrix for A.
It is real-valued for the real Schur decomposition.
sdim : int
If and only if sorting was requested, a third return value will
contain the number of eigenvalues satisfying the sort condition.
Raises
------
LinAlgError
Error raised under three conditions:
1. The algorithm failed due to a failure of the QR algorithm to
compute all eigenvalues
2. If eigenvalue sorting was requested, the eigenvalues could not be
reordered due to a failure to separate eigenvalues, usually because
of poor conditioning
3. If eigenvalue sorting was requested, roundoff errors caused the
leading eigenvalues to no longer satisfy the sorting condition
See also
--------
rsf2csf : Convert real Schur form to complex Schur form
signm(a, disp=True)
Matrix sign function.
Extension of the scalar sign(x) to matrices.
Parameters
----------
A : array, shape(M,M)
Matrix at which to evaluate the sign function
disp : boolean
Print warning if error in the result is estimated large
instead of returning estimated error. (Default: True)
Returns
-------
sgnA : array, shape(M,M)
Value of the sign function at A
(if disp == False)
errest : float
1-norm of the estimated error, ||err||_1 / ||A||_1
Examples
--------
>>> from scipy.linalg import signm, eigvals
>>> a = [[1,2,3], [1,2,1], [1,1,1]]
>>> eigvals(a)
array([ 4.12488542+0.j, -0.76155718+0.j, 0.63667176+0.j])
>>> eigvals(signm(a))
array([-1.+0.j, 1.+0.j, 1.+0.j])
sinhm(A)
Compute the hyperbolic matrix sine.
This routine uses expm to compute the matrix exponentials.
Parameters
----------
A : array, shape(M,M)
Returns
-------
sinhA : array, shape(M,M)
Hyperbolic matrix sine of A
sinm(A)
Compute the matrix sine.
This routine uses expm to compute the matrix exponentials.
Parameters
----------
A : array, shape(M,M)
Returns
-------
sinA : array, shape(M,M)
Matrix cosine of A
solve(a, b, sym_pos=False, lower=False, overwrite_a=False, overwrite_b=False, debug=False)
Solve the equation ``a x = b`` for ``x``.
Parameters
----------
a : array_like, shape (M, M)
A square matrix.
b : array_like, shape (M,) or (M, N)
Right-hand side matrix in ``a x = b``.
sym_pos : bool
Assume `a` is symmetric and positive definite.
lower : boolean
Use only data contained in the lower triangle of `a`, if `sym_pos` is
true. Default is to use upper triangle.
overwrite_a : bool
Allow overwriting data in `a` (may enhance performance).
Default is False.
overwrite_b : bool
Allow overwriting data in `b` (may enhance performance).
Default is False.
Returns
-------
x : array, shape (M,) or (M, N) depending on `b`
Solution to the system ``a x = b``.
Raises
------
LinAlgError
If `a` is singular.
Examples
--------
Given `a` and `b`, solve for `x`:
>>> a = np.array([[3,2,0],[1,-1,0],[0,5,1]])
>>> b = np.array([2,4,-1])
>>> x = linalg.solve(a,b)
>>> x
array([ 2., -2., 9.])
>>> np.dot(a, x) == b
array([ True, True, True], dtype=bool)
solve_banded((l, u), ab, b, overwrite_ab=False, overwrite_b=False, debug=False)
Solve the equation a x = b for x, assuming a is banded matrix.
The matrix a is stored in ab using the matrix diagonal ordered form::
ab[u + i - j, j] == a[i,j]
Example of ab (shape of a is (6,6), u=1, l=2)::
* a01 a12 a23 a34 a45
a00 a11 a22 a33 a44 a55
a10 a21 a32 a43 a54 *
a20 a31 a42 a53 * *
Parameters
----------
(l, u) : (integer, integer)
Number of non-zero lower and upper diagonals
ab : array, shape (l+u+1, M)
Banded matrix
b : array, shape (M,) or (M, K)
Right-hand side
overwrite_ab : boolean
Discard data in ab (may enhance performance)
overwrite_b : boolean
Discard data in b (may enhance performance)
Returns
-------
x : array, shape (M,) or (M, K)
The solution to the system a x = b
solve_continuous_are(a, b, q, r)
Solves the continuous algebraic Riccati equation, or CARE, defined
as (A'X + XA - XBR^-1B'X+Q=0) directly using a Schur decomposition
method.
Parameters
----------
a : array_like
m x m square matrix
b : array_like
m x n matrix
q : array_like
m x m square matrix
r : array_like
Non-singular n x n square matrix
Returns
-------
x : array_like
Solution (m x m) to the continuous algebraic Riccati equation
Notes
-----
Method taken from:
Laub, "A Schur Method for Solving Algebraic Riccati Equations."
U.S. Energy Research and Development Agency under contract
ERDA-E(49-18)-2087.
http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf
See Also
--------
solve_discrete_are : Solves the discrete algebraic Riccati equation
solve_discrete_are(a, b, q, r)
Solves the disctrete algebraic Riccati equation, or DARE, defined as
(X = A'XA-(A'XB)(R+B'XB)^-1(B'XA)+Q), directly using a Schur decomposition
method.
Parameters
----------
a : array_like
Non-singular m x m square matrix
b : array_like
m x n matrix
q : array_like
m x m square matrix
r : array_like
Non-singular n x n square matrix
Returns
-------
x : array_like
Solution to the continuous Lyapunov equation
Notes
-----
Method taken from:
Laub, "A Schur Method for Solving Algebraic Riccati Equations."
U.S. Energy Research and Development Agency under contract
ERDA-E(49-18)-2087.
http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf
See Also
--------
solve_continuous_are : Solves the continuous algebraic Riccati equation
solve_discrete_lyapunov(a, q)
Solves the Discrete Lyapunov Equation (A'XA-X=-Q) directly.
Parameters
----------
a : array_like
A square matrix
q : array_like
Right-hand side square matrix
Returns
-------
x : array_like
Solution to the continuous Lyapunov equation
Notes
-----
Algorithm is based on a direct analytical solution from:
Hamilton, James D. Time Series Analysis, Princeton: Princeton University
Press, 1994. 265. Print.
http://www.scribd.com/doc/20577138/Hamilton-1994-Time-Series-Analysis
solve_lyapunov(a, q)
Solves the continuous Lyapunov equation (AX + XA^H = Q) given the values
of A and Q using the Bartels-Stewart algorithm.
Parameters
----------
a : array_like
A square matrix
q : array_like
Right-hand side square matrix
Returns
-------
x : array_like
Solution to the continuous Lyapunov equation
Notes
-----
Because the continuous Lyapunov equation is just a special form of the
Sylvester equation, this solver relies entirely on solve_sylvester for a
solution.
See Also
--------
solve_sylvester : computes the solution to the Sylvester equation
solve_sylvester(a, b, q)
Computes a solution (X) to the Sylvester equation (AX + XB = Q).
Parameters
----------
a : array, shape (M, M)
Leading matrix of the Sylvester equation
b : array, shape (N, N)
Trailing matrix of the Sylvester equation
q : array, shape (M, N)
Right-hand side
Returns
-------
x : array, shape (M, N)
The solution to the Sylvester equation.
Raises
------
LinAlgError
If solution was not found
Notes
-----
Computes a solution to the Sylvester matrix equation via the Bartels-
Stewart algorithm. The A and B matrices first undergo Schur
decompositions. The resulting matrices are used to construct an
alternative Sylvester equation (``RY + YS^T = F``) where the R and S
matrices are in quasi-triangular form (or, when R, S or F are complex,
triangular form). The simplified equation is then solved using
``*TRSYL`` from LAPACK directly.
solve_triangular(a, b, trans=0, lower=False, unit_diagonal=False, overwrite_b=False, debug=False)
Solve the equation `a x = b` for `x`, assuming a is a triangular matrix.
Parameters
----------
a : array, shape (M, M)
b : array, shape (M,) or (M, N)
lower : boolean
Use only data contained in the lower triangle of a.
Default is to use upper triangle.
trans : {0, 1, 2, 'N', 'T', 'C'}
Type of system to solve:
======== =========
trans system
======== =========
0 or 'N' a x = b
1 or 'T' a^T x = b
2 or 'C' a^H x = b
======== =========
unit_diagonal : boolean
If True, diagonal elements of A are assumed to be 1 and
will not be referenced.
overwrite_b : boolean
Allow overwriting data in b (may enhance performance)
Returns
-------
x : array, shape (M,) or (M, N) depending on b
Solution to the system a x = b
Raises
------
LinAlgError
If a is singular
Notes
-----
.. versionadded:: 0.9.0
solveh_banded(ab, b, overwrite_ab=False, overwrite_b=False, lower=False)
Solve equation a x = b. a is Hermitian positive-definite banded matrix.
The matrix a is stored in ab either in lower diagonal or upper
diagonal ordered form:
ab[u + i - j, j] == a[i,j] (if upper form; i <= j)
ab[ i - j, j] == a[i,j] (if lower form; i >= j)
Example of ab (shape of a is (6,6), u=2)::
upper form:
* * a02 a13 a24 a35
* a01 a12 a23 a34 a45
a00 a11 a22 a33 a44 a55
lower form:
a00 a11 a22 a33 a44 a55
a10 a21 a32 a43 a54 *
a20 a31 a42 a53 * *
Cells marked with * are not used.
Parameters
----------
ab : array, shape (u + 1, M)
Banded matrix
b : array, shape (M,) or (M, K)
Right-hand side
overwrite_ab : boolean
Discard data in ab (may enhance performance)
overwrite_b : boolean
Discard data in b (may enhance performance)
lower : boolean
Is the matrix in the lower form. (Default is upper form)
Returns
-------
x : array, shape (M,) or (M, K)
The solution to the system a x = b
sqrtm(A, disp=True)
Matrix square root.
Parameters
----------
A : array, shape(M,M)
Matrix whose square root to evaluate
disp : boolean
Print warning if error in the result is estimated large
instead of returning estimated error. (Default: True)
Returns
-------
sgnA : array, shape(M,M)
Value of the sign function at A
(if disp == False)
errest : float
Frobenius norm of the estimated error, ||err||_F / ||A||_F
Notes
-----
Uses algorithm by Nicholas J. Higham
svd(a, full_matrices=True, compute_uv=True, overwrite_a=False)
Singular Value Decomposition.
Factorizes the matrix a into two unitary matrices U and Vh, and
a 1-D array s of singular values (real, non-negative) such that
``a == U*S*Vh``, where S is a suitably shaped matrix of zeros with
main diagonal s.
Parameters
----------
a : ndarray
Matrix to decompose, of shape ``(M,N)``.
full_matrices : bool, optional
If True, `U` and `Vh` are of shape ``(M,M)``, ``(N,N)``.
If False, the shapes are ``(M,K)`` and ``(K,N)``, where
``K = min(M,N)``.
compute_uv : bool, optional
Whether to compute also `U` and `Vh` in addition to `s`.
Default is True.
overwrite_a : bool, optional
Whether to overwrite `a`; may improve performance.
Default is False.
Returns
-------
U : ndarray
Unitary matrix having left singular vectors as columns.
Of shape ``(M,M)`` or ``(M,K)``, depending on `full_matrices`.
s : ndarray
The singular values, sorted in non-increasing order.
Of shape (K,), with ``K = min(M, N)``.
Vh : ndarray
Unitary matrix having right singular vectors as rows.
Of shape ``(N,N)`` or ``(K,N)`` depending on `full_matrices`.
For ``compute_uv = False``, only `s` is returned.
Raises
------
LinAlgError
If SVD computation does not converge.
See also
--------
svdvals : Compute singular values of a matrix.
diagsvd : Construct the Sigma matrix, given the vector s.
Examples
--------
>>> from scipy import linalg
>>> a = np.random.randn(9, 6) + 1.j*np.random.randn(9, 6)
>>> U, s, Vh = linalg.svd(a)
>>> U.shape, Vh.shape, s.shape
((9, 9), (6, 6), (6,))
>>> U, s, Vh = linalg.svd(a, full_matrices=False)
>>> U.shape, Vh.shape, s.shape
((9, 6), (6, 6), (6,))
>>> S = linalg.diagsvd(s, 6, 6)
>>> np.allclose(a, np.dot(U, np.dot(S, Vh)))
True
>>> s2 = linalg.svd(a, compute_uv=False)
>>> np.allclose(s, s2)
True
svdvals(a, overwrite_a=False)
Compute singular values of a matrix.
Parameters
----------
a : ndarray
Matrix to decompose, of shape ``(M, N)``.
overwrite_a : bool, optional
Whether to overwrite `a`; may improve performance.
Default is False.
Returns
-------
s : ndarray
The singular values, sorted in decreasing order.
Of shape ``(K,)``, with``K = min(M, N)``.
Raises
------
LinAlgError
If SVD computation does not converge.
See also
--------
svd : Compute the full singular value decomposition of a matrix.
diagsvd : Construct the Sigma matrix, given the vector s.
tanhm(A)
Compute the hyperbolic matrix tangent.
This routine uses expm to compute the matrix exponentials.
Parameters
----------
A : array, shape(M,M)
Returns
-------
tanhA : array, shape(M,M)
Hyperbolic matrix tangent of A
tanm(A)
Compute the matrix tangent.
This routine uses expm to compute the matrix exponentials.
Parameters
----------
A : array, shape(M,M)
Returns
-------
tanA : array, shape(M,M)
Matrix tangent of A
toeplitz(c, r=None)
Construct a Toeplitz matrix.
The Toeplitz matrix has constant diagonals, with c as its first column
and r as its first row. If r is not given, ``r == conjugate(c)`` is
assumed.
Parameters
----------
c : array_like
First column of the matrix. Whatever the actual shape of `c`, it
will be converted to a 1-D array.
r : array_like
First row of the matrix. If None, ``r = conjugate(c)`` is assumed;
in this case, if c[0] is real, the result is a Hermitian matrix.
r[0] is ignored; the first row of the returned matrix is
``[c[0], r[1:]]``. Whatever the actual shape of `r`, it will be
converted to a 1-D array.
Returns
-------
A : array, shape (len(c), len(r))
The Toeplitz matrix. Dtype is the same as ``(c[0] + r[0]).dtype``.
See also
--------
circulant : circulant matrix
hankel : Hankel matrix
Notes
-----
The behavior when `c` or `r` is a scalar, or when `c` is complex and
`r` is None, was changed in version 0.8.0. The behavior in previous
versions was undocumented and is no longer supported.
Examples
--------
>>> from scipy.linalg import toeplitz
>>> toeplitz([1,2,3], [1,4,5,6])
array([[1, 4, 5, 6],
[2, 1, 4, 5],
[3, 2, 1, 4]])
>>> toeplitz([1.0, 2+3j, 4-1j])
array([[ 1.+0.j, 2.-3.j, 4.+1.j],
[ 2.+3.j, 1.+0.j, 2.-3.j],
[ 4.-1.j, 2.+3.j, 1.+0.j]])
tri(N, M=None, k=0, dtype=None)
Construct (N, M) matrix filled with ones at and below the k-th diagonal.
The matrix has A[i,j] == 1 for i <= j + k
Parameters
----------
N : integer
The size of the first dimension of the matrix.
M : integer or None
The size of the second dimension of the matrix. If `M` is None,
`M = N` is assumed.
k : integer
Number of subdiagonal below which matrix is filled with ones.
`k` = 0 is the main diagonal, `k` < 0 subdiagonal and `k` > 0
superdiagonal.
dtype : dtype
Data type of the matrix.
Returns
-------
A : array, shape (N, M)
Examples
--------
>>> from scipy.linalg import tri
>>> tri(3, 5, 2, dtype=int)
array([[1, 1, 1, 0, 0],
[1, 1, 1, 1, 0],
[1, 1, 1, 1, 1]])
>>> tri(3, 5, -1, dtype=int)
array([[0, 0, 0, 0, 0],
[1, 0, 0, 0, 0],
[1, 1, 0, 0, 0]])
tril(m, k=0)
Make a copy of a matrix with elements above the k-th diagonal zeroed.
Parameters
----------
m : array
Matrix whose elements to return
k : integer
Diagonal above which to zero elements.
k == 0 is the main diagonal, k < 0 subdiagonal and k > 0 superdiagonal.
Returns
-------
A : array, shape m.shape, dtype m.dtype
Examples
--------
>>> from scipy.linalg import tril
>>> tril([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], -1)
array([[ 0, 0, 0],
[ 4, 0, 0],
[ 7, 8, 0],
[10, 11, 12]])
triu(m, k=0)
Make a copy of a matrix with elements below the k-th diagonal zeroed.
Parameters
----------
m : array
Matrix whose elements to return
k : integer
Diagonal below which to zero elements.
k == 0 is the main diagonal, k < 0 subdiagonal and k > 0 superdiagonal.
Returns
-------
A : array, shape m.shape, dtype m.dtype
Examples
--------
>>> from scipy.linalg import triu
>>> triu([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], -1)
array([[ 1, 2, 3],
[ 4, 5, 6],
[ 0, 8, 9],
[ 0, 0, 12]])
DATA
__all__ = ['LinAlgError', 'all_mat', 'basic', 'blas', 'block_diag', 'c...
__version__ = '0.4.9'
VERSION
0.4.9
יתן לכם פירוט מלא של הספריה ויכולותיה.
הספריה matplotlib
matplotlib logo
ספריה ליצירת גרפים ותרשימים. בעיקר בדו-ממד.
לטעמי, הדרך הכי טובה ללמוד להשתמש בה הוא לסרוק את גלריית הדוגמאות ולהעתיק את הגרף שנראה הכי כמו-שאתם-צריכים.
In [15]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib
WARNING: pylab import has clobbered these variables: ['draw_if_interactive', 'new_figure_manager'] `%pylab --no-import-all` prevents importing * from pylab and numpy
In [21]:
import numpy as np
from matplotlib.pyplot import *
x = np.linspace(0,10,100)
y = np.sin(x)
plot(x,y,".-",label="sine")
#plot(y**2,"o-",label="sine^2")
legend(loc="best")
#ylim([0,0.5])
xlabel("x")
ylabel("y")
Out[21]:
<matplotlib.text.Text at 0x444ca50>
In [3]:
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0, 1)
y = np.sin(4 * np.pi * x) * np.exp(-5 * x)
plt.fill(x, y, 'r')
plt.grid(True)
plt
Out[3]:
<module 'matplotlib.pyplot' from '/usr/lib/python2.7/dist-packages/matplotlib/pyplot.pyc'>
In [4]:
import matplotlib.pyplot as plt ; import random
fontsizes = [8, 16, 24, 32]
def example_plot(ax):
ax.plot([1, 2])
ax.set_xlabel('x-label', fontsize=random.choice(fontsizes))
ax.set_ylabel('y-label', fontsize=random.choice(fontsizes))
ax.set_title('Title', fontsize=random.choice(fontsizes))
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2, ncols=2)
example_plot(ax1)
example_plot(ax2)
example_plot(ax3)
example_plot(ax4)
plt.tight_layout()
In [5]:
from pylab import *
import matplotlib.cbook as cbook
# data are 256x256 16 bit integers
dfile = cbook.get_sample_data('s1045.ima.gz')
im = np.fromstring(dfile.read(), np.uint16).astype(float)
im.shape = 256, 256
imshow(im, cmap=cm.jet)
colorbar()
show()
In [6]:
from numpy.random import uniform, seed
from matplotlib.mlab import griddata
import matplotlib.pyplot as plt
import numpy as np
# make up data.
seed(0)
npts = 200
x = uniform(-2,2,npts)
y = uniform(-2,2,npts)
z = x*np.exp(-x**2-y**2)
# define grid.
xi = np.linspace(-2.1,2.1,100)
yi = np.linspace(-2.1,2.1,200)
# grid the data.
zi = griddata(x,y,z,xi,yi,interp='linear')
# contour the gridded data, plotting dots at the nonuniform data points.
CS = plt.contour(xi,yi,zi,15,linewidths=0.5,colors='k')
CS = plt.contourf(xi,yi,zi,15,cmap=plt.cm.rainbow,
vmax=abs(zi).max(), vmin=-abs(zi).max())
plt.colorbar() # draw colorbar
# plot data points.
plt.scatter(x,y,marker='o',c='b',s=5,zorder=10)
plt.xlim(-2,2)
plt.ylim(-2,2)
plt.title('griddata test (%d points)' % npts)
plt.show()
In [10]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def lorenz(x, y, z, s=10, r=28, b=2.667) :
x_dot = s*(y - x)
y_dot = r*x - y - x*z
z_dot = x*y - b*z
return x_dot, y_dot, z_dot
dt = 0.01
stepCnt = 10000
# Need one more for the initial values
xs = np.empty((stepCnt + 1,))
ys = np.empty((stepCnt + 1,))
zs = np.empty((stepCnt + 1,))
# Setting initial values
xs[0], ys[0], zs[0] = (0., 1., 1.05)
# Stepping through "time".
for i in range(stepCnt) :
# Derivatives of the X, Y, Z state
x_dot, y_dot, z_dot = lorenz(xs[i], ys[i], zs[i])
xs[i + 1] = xs[i] + (x_dot * dt)
ys[i + 1] = ys[i] + (y_dot * dt)
zs[i + 1] = zs[i] + (z_dot * dt)
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(xs, ys, zs)
ax.set_xlabel("X Axis")
ax.set_ylabel("Y Axis")
ax.set_zlabel("Z Axis")
ax.set_title("Lorenz Attractor")
plt.show()
matplotlib, xkcd style
In [11]:
import numpy as np
import matplotlib.pyplot as plt
plt.xkcd()
plt.plot(sin(linspace(0,10)),label="line")
plt.title("xkcd style")
plt.legend(loc="best")
plt.xlabel("x axis")
Out[11]:
<matplotlib.text.Text at 0x4a4bd90>
סביבת העבודה, IPython
IPython logo
סביבת עבודה וסט כלים לעבודה אינטראקטיבית. אני לא אדבר עליה באופן פרטני, מאחר ורוב ההרצאה היא פרסומת אחת גדולה ל-IPython
חוץ מסביבת העבודה האינטראקטיבית, IPython כוללת גם כלים למקבול הריצה, אבל לא נתעקב עליהם כאן.
הספריה Pandas
Pandas
הספריה הזו מאפשרת טיפול בנתונים כטבלה: עם עמודות מובחנות ושורות מובחנות וביצוע פעולות "חכמות" עליהם.
פעולות מהירות ומתוחכמות על טבלאות שוכנות זכרון.
In [60]:
import pandas as pd
df = pd.read_csv("/home/ronen/Dropbox/master/alab/squid_data/R01C100_01_ovend.res", sep="\t", header= 1 )
df[:5]
Out[60]:
| time[s] | Tb[K]before | field[T] | average Z-Moment | Ta[K] | Tb[K]after | STD | Z-Scan1 | Z-Scan2 | Z-Scan3 | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 3434870100 | 8.7873 | 0.0002 | -0 | 6.0252 | 8.5841 | 1.773417e-10 | -0 | -0 | -0 |
| 1 | 3434870200 | 8.3708 | 0.0002 | -0 | 5.9130 | 8.1694 | 5.117923e-10 | -0 | -0 | -0 |
| 2 | 3434870300 | 8.1145 | 0.0002 | -0 | 5.8433 | 7.9177 | 2.151293e-10 | -0 | -0 | -0 |
| 3 | 3434870300 | 7.8748 | 0.0002 | -0 | 6.2381 | 7.8609 | 2.835875e-10 | -0 | -0 | -0 |
| 4 | 3434870300 | 7.9014 | 0.0002 | -0 | 6.6219 | 8.1232 | 2.930437e-10 | -0 | -0 | -0 |
In [61]:
df.plot(x="Tb[K]before", y="average Z-Moment", style="*")
Out[61]:
<matplotlib.axes.AxesSubplot at 0x4426650>
הספריה Sympy
sympy logo
מספר דוגמאות לשימוש ב-Sympy
אלגברה בסיסית
אתחול: הרבה פעמים מבוצע בצורה אוטומטית.
In [44]:
from IPython.display import display
from sympy.interactive import printing
printing.init_printing()
from __future__ import division
import sympy as sym
from sympy import *
x, y, z = symbols("x y z")
k, m, n = symbols("k m n", integer=True)
f, g, h = map(Function, 'fgh')
In [45]:
Rational(3,2)*pi + exp(I*x) / (x**2 + y)
Out[45]:
$$\frac{3 \pi}{2} + \frac{e^{i x}}{x^{2} + y}$$
In [46]:
eq = ((x+y)**2 * (x+1))
eq
Out[46]:
$$\left(x + 1\right) \left(x + y\right)^{2}$$
In [47]:
expand(eq)
Out[47]:
$$x^{3} + 2 x^{2} y + x^{2} + x y^{2} + 2 x y + y^{2}$$
In [25]:
from sympy import *
eq = Eq(x**3 + 2*x**2 + 4*x + 8, 0)
eq
--------------------------------------------------------------------------- SympifyError Traceback (most recent call last) <ipython-input-25-113a14de44d4> in <module>() 1 from sympy import * ----> 2 eq = Eq(x**3 + 2*x**2 + 4*x + 8, 0) 3 eq 4 5 sympy.latex(eq) /usr/local/lib/python2.7/dist-packages/sympy/core/relational.pyc in Eq(a, b) /usr/local/lib/python2.7/dist-packages/sympy/core/relational.pyc in __new__(cls, lhs, rhs, rop, **assumptions) /usr/local/lib/python2.7/dist-packages/sympy/core/sympify.pyc in _sympify(a) /usr/local/lib/python2.7/dist-packages/sympy/core/sympify.pyc in sympify(a, locals, convert_xor, strict, rational) SympifyError: SympifyError: array([ 8. , 8.4254771 , 8.89795001, 9.42360242, 10.00861796, 10.65918031, 11.38147313, 12.18168008, 13.06598481, 14.040571 , 15.11162229, 16.28532236, 17.56785486, 18.96540345, 20.48415179, 22.13028355, 23.90998239, 25.82943196, 27.89481593, 30.11231796, 32.4881217 , 35.02841083, 37.739369 , 40.62717987, 43.6980271 , 46.95809436, 50.4135653 , 54.07062359, 57.93545289, 62.01423685, 66.31315914, 70.83840342, 75.59615335, 80.59259259, 85.83390481, 91.32627365, 97.07588279, 103.08891589, 109.3715566 , 115.92998859, 122.77039552, 129.89896104, 137.32186883, 145.04530253, 153.07544582, 161.41848234, 170.08059578, 179.06796977, 188.38678799, 198.0432341 , 208.04349175, 218.39374461, 229.10017634, 240.1689706 , 251.60631104, 263.41838134, 275.61136516, 288.19144614, 301.16480796, 314.53763428, 328.31610875, 342.50641503, 357.1147368 , 372.1472577 , 387.6101614 , 403.50963157, 419.85185185, 436.64300592, 453.88927743, 471.59685004, 489.77190742, 508.42063323, 527.54921112, 547.16382476, 567.27065781, 587.87589393, 608.98571677, 630.60631001, 652.74385731, 675.40454231, 698.59454869, 722.32006011, 746.58726022, 771.40233268, 796.77146117, 822.70082933, 849.19662084, 876.26501934, 903.9122085 , 932.14437199, 960.96769346, 990.38835658, 1020.412545 , 1051.04644239, 1082.2962324 , 1114.1680987 , 1146.66822495, 1179.80279481, 1213.57799194, 1248. ])
In [49]:
solve(eq, x)
Out[49]:
$$\begin{bmatrix}-2, & - 2 i, & 2 i\end{bmatrix}$$
In [50]:
myFunc = cos(x**2)**2 / (1+x); myFunc
Out[50]:
$$\frac{\cos^{2}{\left (x^{2} \right )}}{x + 1}$$
In [51]:
diff(myFunc, x)
Out[51]:
$$- \frac{4 x \cos{\left (x^{2} \right )}}{x + 1} \sin{\left (x^{2} \right )} - \frac{\cos^{2}{\left (x^{2} \right )}}{\left(x + 1\right)^{2}}$$
In [52]:
eqn = Eq(Derivative(f(x),x,x) + 9*f(x), 1)
display(eqn)
dsolve(eqn, f(x))
$$9 f{\left (x \right )} + \frac{d^{2}}{d x^{2}} f{\left (x \right )} = 1$$
Out[52]:
$$f{\left (x \right )} = C_{1} \sin{\left (3 x \right )} + C_{2} \cos{\left (3 x \right )} + \frac{1}{9}$$
פונקציות מיוחדות
אפשר לשחק גם עם דברים יותר קוונטיים, למשל, פונקציות גל של אטום מימן.
In [58]:
from sympy import Eq, Integral, oo, pprint, symbols
from sympy.physics.hydrogen import R_nl
print("Hydrogen radial wavefunctions:")
a, r = symbols("a r")
print("R_{21}:")
R_nl(2, 1, a, r)
Hydrogen radial wavefunctions:
R_{21}:
Out[58]:
$$\frac{a r}{12} \sqrt{6} \sqrt{r^{3}} e^{- \frac{a r}{2}}$$
In [59]:
print("Normalization:")
i = Integral(R_nl(2, 1, 1, r)**2 * r**2, (r, 0, oo))
Eq(i, i.doit())
Normalization:
Out[59]:
$$\int_{0}^{\infty} \frac{r^{7}}{24} e^{- r}\, dr = 210$$
ממש קוונטים - תנע זוויתי
In [54]:
from sympy.physics.quantum.spin import Jminus, Jx, Jz, J2, J2Op, JzKet, JzKetCoupled, Rotation, WignerD, couple, uncouple
from sympy.physics.quantum import Dagger, hbar, qapply, represent, TensorProduct
from sympy import factor, pi, S, Sum, symbols
In [55]:
jz = JzKet(1,1); jz
Out[55]:
$${\left|1,1\right\rangle }$$
In [56]:
represent(jz)
Out[56]:
$$\left[\begin{matrix}1\\0\\0\end{matrix}\right]$$
In [24]:
ip = Dagger(jz) * jz; ip
Out[24]:
$$\left\langle 1,1 \right. {\left|1,1\right\rangle }$$
In [25]:
ip.doit()
Out[25]:
$$1$$
In [26]:
Jz * jz
Out[26]:
$$J_z {\left|1,1\right\rangle }$$
In [27]:
qapply(Jz * jz)
Out[27]:
$$\hbar {\left|1,1\right\rangle }$$
In [28]:
Jminus * jz
Out[28]:
$$J_- {\left|1,1\right\rangle }$$
In [29]:
qapply(Jminus * jz)
Out[29]:
$$\sqrt{2} \hbar {\left|1,0\right\rangle }$$
In [30]:
j, m = symbols('j m')
jz = JzKet(j, m); jz
Out[30]:
$${\left|j,m\right\rangle }$$
In [31]:
J2 * jz
Out[31]:
$$J^2 {\left|j,m\right\rangle }$$
In [32]:
qapply(J2 * jz)
Out[32]:
$$\hbar^{2} j^{2} {\left|j,m\right\rangle } + \hbar^{2} j {\left|j,m\right\rangle }$$