SciPy Linalg

SciPy is built using the optimized ATLAS LAPACK and BLAS libraries. It has very fast linear algebra capabilities. All of these linear algebra routines expect an object that can be converted to a two-dimensional array. The output of these routines is also a two-dimensional array.

SciPy.linalg vs. NumPy.linalg

A scipy.linalg includes all the functions in numpy.linalg. In addition, scipy.linalg has some other advanced functions that are not found in numpy.linalg. Another benefit of using scipy.linalg over numpy.linalg is that it is always compiled with BLAS/LAPACK support, whereas this is optional for NumPy. Therefore, the SciPy version may be faster, depending on how NumPy is installed.

Linear Equations

scipy.linalg.solve solves the linear equation a * x + b * y = Z for unknown x and y values.

As an example, suppose we need to solve the following homology equation.

x + 3y + 5z = 10

2x + 5y + z = 8

2x + 3y + 8z = 3begin{bmatrix} 1 & 3 & 5 2 & 5 & 1 2 & 3 & 8end{bmatrix}^{-1} begin{bmatrix} 10 8 3end{bmatrix} = frac{1}{25}.begin{bmatrix} -232 129 19 end{bmatrix} =begin{bmatrix} -9.28 5.16 0.76 `end{bmatrix}.

However, it is better to use the linalg.solve command, which is faster and more stable for numerical computations.

The solve function takes two inputs, ‘a’ and ‘b’, where ‘a’ represents the coefficients and ‘b’ represents the respective right-hand side values, and returns an array of the solution.

Let’s consider the following example.

#importing the scipy and numpy packages
from scipy import linalg
import numpy as np

#Declaring the numpy arrays
a = np.array([[3, 2, 0], [1, -1, 0], [0, 5, 1]])
b = np.array([2, 4, -1])

#Passing the values to the solve function
x = linalg.solve(a, b)

#printing the result array
print x

The above program will produce the following output.

array([ 2., -2., 9.])

Finding the Determinant

The determinant of a square matrix A, usually denoted |A|, is a common quantity in linear algebra. In SciPy, this is calculated using the det() function. It takes a matrix as input and returns a scalar value.

Let’s take a look at the following example.

#Importing the scipy and numpy packages
from scipy import linalg
import numpy as np

#Declaring the numpy array
A = np.array([[1,2],[3,4]])

#Passing the values to the det function
x = linalg.det(A)

#Printing the result
print x

The above program will produce the following output.

-2.0

Eigenvalues and Eigenvectors

The eigenvalue-eigenvector problem is one of the most commonly used linear algebra operations. We can find the eigenvalues (λ) and corresponding eigenvectors (v) of a square matrix (A) by considering the following relationship.

Av = λv

scipy.linalg.eig computes eigenvalues from a general or generalized eigenvalue problem. This function returns both the eigenvalue and the eigenvector.

Let’s consider the following example.

#importing the scipy and numpy packages
from scipy import linalg
import numpy as np

#Declaring the numpy array
A = np.array([[1,2],[3,4]])

#Passing the values to the eig function
l, v = linalg.eig(A)

#printing the result for eigen values
print l

#printing the result for eigen vectors
print v

The above program will produce the following output.

array([-0.37228132+0.j, 5.37228132+0.j]) #--Eigen Values
array([[-0.82456484, -0.41597356], #--Eigen Vectors
[ 0.56576746, -0.90937671]])

Singular Value Decomposition

Singular Value Decomposition (SVD) can be considered an extension of the eigenvalue problem for non-square matrices.

scipy.linalg.svd decomposes a matrix ‘a’ into two unit matrices ‘U’ and ‘Vh’ and a one-dimensional array ‘s’ of real, non-negative singular values, such that a == USVh, where ‘S’ is a zero matrix of appropriate shape with ‘s’ as the main diagonal.

Let’s consider the following example.

#importing the scipy and numpy packages
from scipy import linalg
import numpy as np

#Declaring the numpy array
a = np.random.randn(3, 2) + 1.j*np.random.randn(3, 2)

#Passing the values to the eig function
U, s, Vh = linalg.svd(a)

# printing the result
print U, Vh, s

The above program will produce the following output.

(
   array([
      [0.54828424-0.23329795j, -0.38465728+0.01566714j,
      -0.18764355+0.67936712j],
      [-0.27123194-0.5327436j , -0.57080163-0.00266155j,
      -0.39868941-0.39729416j],
      [ 0.34443818+0.4110186j , -0.47972716+0.54390586j,      0.25028608-0.35186815j]
   ]),

   array([ 3.25745379, 1.16150607]),

   array([
      [-0.35312444+0.j , 0.32400401+0.87768134j],
      [-0.93557636+0.j , -0.12229224-0.33127251j]
   ])
)