The QR decomposition (also known as QR factorization) is the decomposition of a matrix A into a product of an orthogonal matrix Q and upper triangular matrix R. A=QR and Q^T Q=I [A:4].

The decomposition is unique if A is a real, square, and invertible matrix. In the case of a rectangle matrix A, m by n with m > n the decomposition is implemented as the dot product of two vector of matrices: A = [Q₁ , Q₂ ].[R₁ , R₂ ]^T where Q₁ is an m by n matrix, Q₂ is an m by n matrix, R₁ is n by n and an upper triangle matrix, R₂ is an m by n null matrix.

QR decomposition is a reliable method for solving large systems of linear equations for which the number of equations (rows) exceeds the number of variables (columns). Its asymptotic computational time complexity for a training set of m dimension and n observations is O(mn² -n³ /3).

It is used to minimize the loss function for ordinary least squares regression, which is covered in the Ordinary least squares regression section in Chapter 9, Regression and Regularization.

LU factorization is a technique to solve a matrix equation A.x = b, where A is a nonsingular matrix and x and b are two vectors. The technique consists of decomposing the original matrix A as the product of simple matrix A= A₁ A₂ …A_n.

LDL decomposition for real matrices defines a real positive matrix A as the product of a lower unit triangular matrix L, a diagonal matrix D, and the transposed matrix of L, L^T A=LDL^T.

The Cholesky factorization (also known as Cholesky decomposition) of real matrices is a special case of LU factorization [A:4]. It decomposes a positive-definite matrix A into a product of lower triangle matrix L and its conjugate transpose L^T, A=LL^T.

The asymptotic computational time complexity for the Cholesky factorization is O(mn2), where m is the number of features (model parameters) and n the number of observations. The Cholesky factorization is used in linear least squares, Kalman filter (in the Kalman filter: Recursive algorithm section in Chapter 3, Data Preprocessing), and nonlinear Quasi-Newton optimizer.

The SVD of real matrices defines a m x n real matrix A as the product of a m square real unitary matrix, U a m x n rectangular diagonal matrix S and the transpose V^T matrix of a real matrix, A=USV^T.

The columns of the matrix U and V are the orthogonal bases and the value of the diagonal matrix S is the singular values [A:4]. The asymptotic computational time complexity for the singular value decomposition for n observations and m features is O(mn² -n³ ). The singular value decomposition is used in minimizing the total least squares and solving homogeneous linear equations.

The Eigen-decomposition of a real square matrix A is the canonical factorization as Ax = ?x.

? is the eigenvalue (scalar) corresponding to vector x. The n by n matrix A is then defined as A = QDQ^T . Q is the square matrix that contains the eigenvectors and D is the diagonal matrix wherein the elements are the eigenvalues associated to the eigenvectors [A:5], [A:6]. The Eigen-decomposition is used in Principal Components Analysis (in the Principal Components Analysis section in Chapter 5, Dimension Reduction.)

There are many more Open Source algebraic libraries available to developers as API besides the Apache Commons Math used in Chapter 3, Data Pre-processing, Chapter 5, Dimensional Reduction, and Chapter 6, Naive Bayes Classifiers, and Apache Spark/MLlib in Chapter 17, Apache Spark MLlib.

The Apache Spark/MLlib framework bundles jBlas, Colt, and Breeze. The Iitb framework for conditional random fields uses colt linear algebra components.

There are many optimization techniques that either rely on exclusively on the first order derivative or provide a linear algebra alternative to the computation of the gradient.

Steepest descent

The steepest descent (also known as gradient descent) method is one of the simplest techniques for finding a local minimum of any continuous, differentiable function F or the global minimum for any define, differentiable, convex function [A:9]. The value of a vector or data point x_t+1 at iteration t+1 is computed from the previous value x_t using the gradient F of function F and the slope ?:

The steepest gradient algorithm is used for solving a system of nonlinear equations, minimization of the loss function in the logistic regression (refer to the Numerical optimization section in Chapter 9, Regression and Regularization), support vector classifier (refer to the Non-separable case section in Chapter 12, Kernel Models and Support Vector Machine), and multilayer perceptron (refer to the Training strategy and classification section in Chapter 10, Multilayer Perceptron).

Conjugate gradient

The conjugate gradient solves unconstrained optimization problems and systems of linear equations. It is an alternative to the LU factorization for positive, definite symmetric square matrices. The solution x* of the equation A.x = b is expanded as the weighted summation of n basis orthogonal directions p_i (also known as conjugate directions):

The solution x* is extracted by computing the i^th conjugate vector p_i and then computing the coefficients a_i.

Stochastic gradient descent

The stochastic gradient method is a variant of the steepest descent that minimizes the convex function by defining the objective function F as the sum of differentiable, basis function f_i as follows:

The solution x_t+1 at iteration t+1 is computed from the value x_t at iteration t, the step size (or learning rate) a and the sum of the gradient of the basic functions [A:10]. The stochastic gradient descent is usually faster than other gradient descents or quasi-new methods in converging toward a solution for a convex function. The stochastic gradient descent is used in the logistic regression, support vector machines, and back-propagation neural networks.

Stochastic gradient is particularly suitable for discriminative models with large datasets [A:11]. Spark/MLlib makes extensive use of the stochastic gradient method.

Note

Batch gradient descent

The batch gradient descent is introduced and implemented in the Step 5: Implementing the classifier section under Let's kick the tires in Chapter 1, Getting Started.

Quasi-Newton algorithms are variations on Newton's methods of finding the value of a vector or data point that maximizes or minimizes a function F (1st order derivative is null) [A:12].

Newton's method is a well-known and simple optimization method for finding the solution of equations F(x) = 0 for which F is continuous and 2nd order differentiable. It relies on the Taylor series expansion to approximate the function F with a quadratic approximation on variable ?x = x_t+1 - x_t to compute the value at the next iteration using the first order F' and second order F" derivatives:

Contrary to the Newton method, Quasi-Newton methods do not require that the 2nd order derivative, Hessian matrix, of the objective function be computed. It just has to be approximated [A:13]. There are several approaches to approximate the computation of the Hessian matrix.

BFGS

The Broyden-Fletcher-Goldfarb-Shanno (BFGS) is a Quasi-Newton iterative numerical method to solve unconstrained nonlinear problems. The Hessian matrix H_t+1 at iteration t is approximated using the value of the previous iteration t as H_t+1 =H_t + U_t + V_t applied to the Newton equation for the direction p_t:

The BFGS is used in the minimization of the cost function for the conditional random field and L₁ and L₂ regression.

L-BFGS

The performance of the BFGS algorithm is related to the caching of the approximation of the Hessian matrix in memory (U, V) at the cost of high memory consumption.

The Limited memory BFGS algorithm or L-BFGS is a variant of BFGS that uses a minimum amount of computer RAM. The algorithm maintains the last m incremental updates of the values ?x_t and gradient ?G_t at iteration t, then computes these values for the next step t+1:

It is supported in Apache Commons Math 3.3+, Apache Spark/MLlib 1.0+, Colt 1.0+, and Iiitb CRF libraries. L-BFGS is used in the minimization of the loss function in the Conditional Random Fields (refer to the Conditional random fields section in Chapter 7, Sequential Data Models).

Let's consider the classic minimization of the least squares of a nonlinear function y = F(x, w) with w_i parameters for observations {y, x_i }. The objective is to minimize the sum of the squares of residuals r_i:

Gauss-Newton

The Gauss-Newton technique is a generalization of Newton's method. The technique solves nonlinear least squares by updating the parameters w_t+1 at iteration t+1 using the first order derivative (also known as Jacobian):

The Gauss-Newton algorithm is used in Logistic Regression (refer to the Logistic regression section in Chapter 9, Regression and Regularization).

Levenberg-Marquardt

The Levenberg-Marquardt algorithm is an alternative to the Gauss-Newton for solving nonlinear least squares and curve fitting problems. The method consists of adding the gradient (Jacobian) terms to the residuals r_i to approximate the least squares error:

The algorithm is used in the training of the logistic regression (refer to the Logistic regression section in Chapter 9, Regression and Regularization).

The Lagrange multipliers methodology is an optimization technique to find local optima of a multivariate function subject to equality constraints [A:14]. The problem is stated as follows:

maximize f(x) subject to g(x) = c, c is a constant, x is a variable or features vector.

The methodology introduces a new variable ? to integrate the constraint g into a function, known as the Lagrange function ℒ(x, ?). Let's note the ℒ the gradient of ℒ over the variables x_i and ?. The Lagrange multipliers are computing by maximizing ℒ:

Consider the following example:

Lagrange multipliers are used in minimizing the loss function in the nonseparable case of linear support vector machines (refer to the Non-separable case section of Chapter 12, Kernel Models and Support Vector Machines).