SciPy - lstsq() Function

---

---

SciPy has a function lstsq() to solve linear matrix equations of the form Ax=b with the least squares method. It is applied for over-determined (more equations than unknowns) or under-determined systems and provides a solution x minimizing the error between Ax and b.

The lstsq() method produces key outputs, including residuals, rank, and singular values. Residuals show the gap between the target vector b and the approximation Ax indicating how well the solution matches the data.

The rank of matrix A helps to determine if the matrix has full rank meaning it lacks duplicate rows or columns, which ensures the system's stability.

LAPACK (Linear Algebra PACKage) is a parameter that relates to performing linear algebra operations. These operations include matrix factorizations eigenvalue problems, and solutions to linear systems. This method points out which algorithm solves the least square problem.

Common LAPACK drivers comprise dgesv to solve linear equations dgelss for least-squares problems dgeqrf for QR factorizations, and gelsd as a more reliable method for least-squares problems using a split-and-conquer approach.

Syntax

The syntax for the SciPy lstsq() method is as follows −

.lstsq(a, b, cond=None, overwrite_a=False, overwrite_b=False, check_finite=True, lapack_driver='gelsd')

Parameters

This method accepts the following parameters −

• a (array_like m,n)− The coefficient matrix A.

• b (array_like m or m,p)− The target values B. If B has multiple columns, the solution is computed for each column.

• cond (float, optional) − Threshold for singular values. Singular values smaller than cond are treated as zero.

• overwrite_a (bool, optional) − Allow overwriting data in a to reduce memory usage.

• overwrite_b (bool, optional) − Allow overwriting data in b to reduce memory usage.

• check_finite (bool, optional) − Check whether the input contains only finite numbers. Disabling this improves performance but requires user validation.

• lapack_driver (str, optional) − Specifies the LAPACK driver to use ('gelsd', 'gelss', 'gelsy', or 'gelss'). gelsd is the default driver and generally the most robust.

Return Value

This method returns the following −

• x (n,) or (n, p) ndarray − The least-squares solution to the system Ax=B.

• residuals (p) ndarray or float − The sum of squared residuals for each solution. If the system is underdetermined or the solution is exact, this will be empty.

• rank int − The rank of matrix A, representing the number of linearly independent rows or columns.

• s (min(m, n),) ndarray or NoneThe singular values of A, useful for analyzing the stability and conditioning of the solution.

Example 1: Solving an Overdetermined System

In the below example we have a overdetermined matrix A 32, which has more equations and few unknowns, meaning more rows few columns. For overdetermined matrix there will be no exact solution, so the least square method will find the best solution that minimizes the residuals.

In the below code we have defined coefficient matrix(A) and target matrix (B) and we will find the least-square solution (x) from the equation Ax=B using lstsq() method.

import numpy as np
import scipy.linalg

A = np.array([[1, 1], [1, 2], [1, 3]])
B = np.array([1, 2, 2])

x, residuals, rank, s = scipy.linalg.lstsq(A, B)

print("Solution:", x)
print("Residuals:", residuals)
print("Rank of A:", rank)
print("Singular values of A:", s)

When we run above program, it produces following result −

Solution: [0.66666667 0.5       ]
Residuals: 0.16666666666666677
Rank of A: 2
Singular values of A: [4.07914333 0.60049122]

Example 2: Solving an Underdetermined System

In the below example we have a underdetermined matrix A 23, which has few equations and more unknowns, meaning more columns few rows. For underdetermined matrix there will be infinite solutions.

The code uses lstsq() method to calculate the best solution that satisfies Ax=B exactly. Since the solution satisfies exactly the residual (|ax-b|)^2 will be zero in this case.

import numpy as np
from scipy.linalg import lstsq

# Define the underdetermined coefficient matrix A and target vector B
A = np.array([[1, 0, 1], [0, 1, 1]])
B = np.array([2, 1])

# Compute the least-squares solution
x, residuals, rank, s = lstsq(A, B)

print("Solution:", x)
print("Residuals:", residuals)
print("Rank of A:", rank)

Following is an output of the above code −

Solution: [1.00000000e+00 2.29326682e-16 1.00000000e+00]
Residuals: []
Rank of A: 2

Example 3: Least Square Method with LAPACK Driver Parameter

The lapack_driver='gelsd' uses divide and conquer method for solving least-square problems this improves numerical stability, especially for ill-conditioned matrices (nearly singular).

In the below code we use lstsq() to compute the least-squares solution to the system Ax=b and we will specify the lapack_driver='gelsd' parameter.

import numpy as np
from scipy.linalg import lstsq

# Define the coefficient matrix A and target vector b
A = np.array([[1, 1], [1, 2], [1, 3]])
b = np.array([1, 2, 2])

# Perform least-squares fitting with a specified LAPACK driver
x, residuals, rank, singular_values = lstsq(A, b, lapack_driver='gelsd')

# Output the results
print("Solution:", x)
print("Residuals:", residuals)

Output of the above code is as follows −

Solution: [0.66666667 0.5       ]
Residuals: 0.16666666666666677
Rank of A: 2
Singular values of A: [4.07914333 0.60049122]

Example 4: Solve Noisy Linear Systems with Least Squares Method.

If noise is present in the data, the lstsq() function computes the best fitting solution that has the minimal residual error for approximating Ax=b as close to possible.

Here, we created an overdetermined system with noisy target values of b. Then, using scipy.linalg.lstsq() we calculated the least square solution, along with singular values of A, rank, and residuals.

import numpy as np
from scipy.linalg import lstsq

# Define an overdetermined system with noise
A = np.array([[1, 1], [1, 2], [1, 3]])
b = np.array([1.1, 1.9, 2.1])  # Noisy target values

# Compute the least-squares solution
x, residuals, rank, s = lstsq(A, b)

print("Solution:", x)
print("Residuals:", residuals)
print("Rank of A:", rank)
print("Singular values of A:", s)

Output of the above code is as follows −

Solution: [0.7 0.5]
Residuals: 0.05999999999999997
Rank of A: 2
Singular values of A: [4.07914333 0.60049122]