In this article, I shall go over the topic of arriving at the Vectorized Gradient-Descent formulae for the Cost function of the for Matrix form of training-data Equations. And along with that the Fundamentals of Calculus (especially Partial Derivative) and Matrix Derivatives necessary to understand the process.

So our target of this article is to understand the full Mathematics and the flow behind arriving at the below formulae, which is the Vectorized Gradient of the training-data Matrix

First a Refresher on basic Matrix Algebra

A matrix A over a field K or, simply, a matrix A (when K is implicit) is a rectangular array of scalars usually presented in the following form:

The rows of such a matrix A are the m horizontal lists of scalars:

and the columns of A are the n vertical lists of scalars:

A matrix with m rows and n columns is called an m by n matrix, written m*n. The pair of numbers m and n is called the size of the matrix. Two matrices A and B are equal, written A = B, if they have the same size and if corresponding elements are equal. Thus, the equality of two m * n matrices is equivalent to a system of mn equalities, one for each corresponding pair of elements.
A matrix with only one row is called a row matrix or row vector, and a matrix with only one column is called a column matrix or column vector. A matrix whose entries are all zero is called a zero matrix and will usually be denoted by 0.

Matrix Multiplication

The below image is taken from Khan Academy’s excellent linear algebra course.

In above, each entry in the product matrix is the dot product of a row in the first matrix and a column in the second matrix

More explanation for higher dimension case — If the product AB = C is defined, where C is denoted by [cij], then the
element cij is obtained by multiplying the elements in the ith row of A by the corresponding elements in the jth column of B and adding. Thus, if A has order k * n, and B has order n * p then

then c11 is obtained by multiplying the elements in the first row of A by the corresponding elements in the first column of B and adding; hence,

The element c12 is found by multiplying the elements in the first row of A by the corresponding elements in the second column of B and adding; hence,

The element ckp below is obtained by multiplying the elements in the kth row of A by the corresponding elements in the pth column of B and adding; hence,

There are four simple rules that will help us in multiplying matrices, listed here

1. Firstly, we can only multiply two matrices when the number of columns in
matrix A is equal to the number of rows in matrix B.

2. Secondly, the first row of matrix A multiplied by the first column of matrix B
gives us the first element in the matrix AB, and so on.

3. Thirdly, when multiplying, order matters — specifically, AB ≠ BA.

4. Lastly, the element at row i, column j is the product of the ith row of matrix A and the jth column of matrix B.

The order in which we multiply matters. We must keep the matrices
in order, but we do have some flexibility. As we can see in the following equation, the parentheses can be moved:

The following three rules also apply for Matrix Operation.

Let's see an example of Matrix multiplication

Another example of Matrix multiplication

Implementing a dot production with numpy

import numpy as np# this will create integer as array-elements
A = np.random.randint(10, size=(2,3))
B = np.random.randint(10, size=(3,2))This will create 2 Matrices as below e.g.
[[9 0 8] [3 5 5]][[5 2]
 [0 7]
 [0 7]
]print(np.dot(A, B))# Output 
[[45 74]
[15 76]]

The Hadamard Product — Other Kinds of Matrix Multiplication

Hadamard multiplication is defined for matrices of the same shape as the multiplication of each element of one matrix by the corresponding element of the other matrix. Hadamard multiplication is often denoted by as below, for two matrices A(n×m) and B(n×m) we have

Refresher on Multivariate Linear Regression

First, start with the Simple Linear Regression (SLR).

Suppose we have 2 equations as below

10 = 2x + 2y
18 = 4x + y

So the Matrix form

So in general Mathematic form for the single independent variable case

So the set of equations for all the observation will be as below

And the Matrix form of which is below

So Y is n * 1 matrix, X is an * 2 matrix, β is 2 * 1 matrix

Multiple Linear Regression (MLR)

Suppose that the response variable Y and at least one predictor variable xi are quantitative. Then the equation for a specific Y value under the MLR model is

for i = 1, . . . , n. Here n is the sample size and the random variable ei is the
ith error. So

Which is in compact notation below

And now in matrix notation, these n sets of equations become

where Y is the vector of the response variable and is an n × 1 vector of dependent variables, X is the matrix of the k independent/explanatory variables (usually the first column is a column of ones for the constant term) and is an n × p matrix of predictors, β is a p × 1 vector of unknown coefficients, and e is an n × 1 vector of unknown errors. Equivalently,

Where
Y : is output vector for n training examples.
X : is matrix of size n*p where each ith row belongs to ith training set.
β : is weight vector of size p for p training features.

Different types of Matrix Differentiation

1-Differentiation with Respect to a Scalar

Differentiation of a structure (vector or matrix, for example) with respect to a scalar is quite simple; it just yields the ordinary derivative of each element of the structure in the same structure. Thus, the derivative of a vector or a matrix with respect to a scalar variable is a vector or a matrix, respectively, of the derivatives of the individual elements.

1.A-Derivatives of Vectors with Respect to Scalars

The derivative of the vector y(x) = (y1 , . . . , yn ) with respect to the scalar x
is the vector

The second or higher derivative of a vector with respect to a scalar is likewise a vector of the derivatives of the individual elements; that is, it is an array of higher rank.

1.B-Derivatives of Matrices with Respect to Scalars

The derivative of the matrix Y(x) defined as below

with respect to the scalar, x is the matrix

The second or higher derivative of a matrix with respect to a scalar is
likewise a matrix of the derivatives of the individual elements.

1.C — Derivatives of Functions with Respect to Scalars

Differentiation of a function of a vector or matrix that is linear in the elements
of the vector or matrix involves just the differentiation of the elements, fol-
lowed by application of the function. For example, the derivative of a trace of
a matrix is just the trace of the derivative of the matrix. On the other hand,
the derivative of the determinant of a matrix is not the determinant of the
derivative of the matrix

1.D — Higher-Order Derivatives with Respect to Scalars

Because differentiation with respect to a scalar does not change the rank of the object (“rank” here means rank of an array or “shape”), higher-order derivatives

with respect to scalars are merely objects of the same rank
whose elements are the higher-order derivatives of the individual elements.

2-Differentiation with Respect to a Vector

Differentiation of a given object with respect to an n-vector yields a vector for each element of the given object. The basic expression for the derivative, from formula

for an arbitrary conformable vector y. The arbitrary y indicates that the derivative is omnidirectional; it is the rate of change of a function of the vector in any direction.

2.A — Derivatives of Scalars with Respect to Vectors; The Gradient

The derivative of a scalar-valued function with respect to a vector is a vector
of the partial derivatives of the function with respect to the elements of the
vector. If f (x) is a scalar function of the vector x = (x1 , . . . , xn ),

if those derivatives exist. This vector is called the gradient of the scalar-valued
function, and is sometimes denoted by ∇f (x)

2.B — Derivatives of Vectors with Respect to Vectors; The Jacobian

The derivative of an m-vector-valued function of an n-vector argument consists of nm scalar derivatives. These derivatives could be put into various structures. Two obvious structures are an n × m matrix and an m × n matrix.

For a function,

we define

to be the n × m matrix, which is the natural extension of ∂/∂x applied to a scalar function. The above the notation is more precise because it indicates that the elements of f correspond to the columns of the result.

if those derivatives exist. This derivative is called the matrix gradient and
is denoted by ∇f for the vector-valued function f . (Note that the ∇
symbol can denote either a vector or a matrix, depending on whether the
function being differentiated is scalar-valued or vector-valued.)

The m × n matrix is called the Jacobian of f and is denoted by Jf as below

General Form of the Jacobian

For arriving at the general Mathematical form of Jacobian I would refer a quite well-recognized Paper in this field.

Let y = f(x) be a vector of m scalar-valued functions that each take a vector x of length n = |x| where |x| is the cardinality (count) of elements in x. Each fi function within f returns a scalar just as in the previous section:

For instance,

So the set of equations are as below

So we have m = n functions and parameters, in this case. Generally speaking, though, the Jacobian matrix is the collection of all m × n possible partial derivatives (m rows and n columns), which is the stack of m gradients with respect to x:

Each of these

is a horizontal n-vector because the partial derivative is with respect to a vector, x, whose length is n = |x|. The width of the Jacobian is n if we’re taking the partial derivative with respect to x because there are n parameters we can wiggle, each potentially changing the function’s value. Therefore, the Jacobian is always m rows for m equations.

Now the Cost Function under Multivariate Linear Regression

The equation for the hypothesis function is as follows

The general notations that I will use for extending the above function

So if we are predicting house-price with the above MLR equation, then θ0 will be the basic/base price of a house, then θ1 as the price per room, θ2 as the price per KM-distance from the nearest Airport.

Using the definition of matrix multiplication, our multivariate hypothesis function can be concisely represented as:

This is a vectorization of our hypothesis function for one training example;

Now, using the fact that for a vector z, we have that

Applying the above identity to the right-hand-side of the Cost function (below)

So now the Cost function takes the following form

And our target is to

Wher the thetas θ are the weights, and the above partial derivative for any weights wj will be as below

So the Gradient-Descent process for Multivariate case becomes

Where θ and x are column vector given by

And that's why we take the transpose of θ to multiply with column-vector x to get the hypothesis (as earlier mentioned in this article)

Refresher — Matrix-Derivative Identities required for the Mathematical Derivation of the Gradient of a Matrix w.r.t. to Vectors

The derivative of a matrix is usually referred to as the gradient and denoted as ∇. Consider a function

That is the Matrix will be as below

Then,

Thus, the gradient ∇Af(A) is itself an m-by-n matrix, whose (i, j)-element is

For example, lets take a look at a very simple case. Suppose

is a 2-by-2 matrix, and the function

is given by

Here, Aij denotes the (i, j) entry of the matrix A. We then have

Matrix Transpose Identities

Each of the below identities can be proved separately mathematically proved.

Another related one, If 𝐴 and 𝐵 are two matrices of the same order, then

The Trace: tr(·) of a Matrix

The sum of the diagonal elements of a square matrix is called the trace of the
matrix. We use the notation “tr(A)” to denote the trace of the matrix A:

Because of the associativity of matrix multiplication, this relation can be
extended as

Now proceed to find the Gradient of the Cost function.

First a Refresher on the Gradient Descent Algorithm

To implement Gradient Descent, you need to compute the gradient of the cost function with regard to each model parameter θj. In other words, you need to calculate how much the cost function will change if you change θj just a little bit. This is called a partial derivative. It is like asking “What is the slope of the mountain under my feet if I face east?”. and then asking the same question facing north.

Now you would recognize the very well-known cost function

And the following the Jacobian identity discussed above the Gradient vector of the cost function will be

Notice that this formula involves calculations over the full training set X, at each Gradient Descent step! This is why the algorithm is called Batch Gradient Descent: it uses the whole batch of training data at every step.

Once you have the gradient vector, which points uphill, just go in the opposite direction to go downhill. This means subtracting ∇θMSE(θ) from θ. This is where the learning rate η comes into play:5 multiply the gradient vector by η to determine the size of the downhill step

Now repeating below section of the Matrix form of the training dataset, from our earlier part of this article —

The general form of multiple linear regression (MLR) model is

for i = 1, . . . , n. Here n is the sample size and the random variable ei is the
ith error. In matrix notation, these n sets of equations become

The Y vector is the response variable and is an n × 1 vector of dependent variables, X is the matrix of the k independent/explanatory variables (usually the first column is a column of ones for the constant term) and is an n × p matrix of predictors, β is a p × 1 vector of unknown coefficients, and e is an n × 1 vector of unknown errors. Equivalently,

Where
Y : is output vector for n training examples.
X : is matrix of size n*p where each ith row belongs to ith training set.
β : is weight vector of size p for p training features.

Note that β in the above is not a scalar, but a vector.

Now we have the RSS defined as

Alternative-1 — Vectorized Calculation of Gradient of our Matrix form training data

Note the above is directly derived from using the identity that for a vector z, we have

Now let,

Then for the following assumption

Therefore,

We have, for each of k=1,…, p

Then for the whole matrix (i.e. the whole set of training data set or the whole set of Hypothesis Equation ), we will get

Which in the final Vector form

Note, that in the last equality, I had to get the Transpose of X because when doing matrix multiplication — that's a dot product of rows of the first matrix to columns of the second matrix. The number of columns of the 1st matrix must equal the number of rows of the 2nd matrix.

So by transposing the p-th column of X ends up being the p-th row of the X-Transposed. Thus, when doing a dot product between the p-th row of X-Transposed with (y — Xβ) it will match perfectly as I am using all of the entries of the p-th column of X

Alternative-2 — And below is an alternative calculation for arriving at the same Vectorized Gradient formulae for training-data in Matrix form.

Here, I am denoting the coefficients with θ or Theta (instead of β that we used above in our Alternative-1 Gradient Calculation — only to make the presentation differentiable)

Again assume we have our Training Set of data as below

Also, let y be the m-dimensional vector containing all the target values from the training set:

And we have the Predicted Value or the Hypothesized value as below

So we can say the below

And now again, we need to use the same vector identity mentioned above, that for a vector z, we have

Using the above we have the below relation for the Cost function

We have already introduced the trace operator of a Matrix, written as “tr.” Now we need to use a couple of more matrix derivatives Identities (that I am just stating below here, and they all have robust Mathematical proofs, the details of which I am not including here).

So below 2 Matrix Derivative Identities hold true and we need to use them to arrive at the Gradient Calculation.

Combining the above two Equations or Identities we derive

So now Final Gradient Calculation will be as below

In the third step above, we used the fact that the trace of a real number is just the real number; the fourth step used the fact that

And the fifth step used below equation that we already mentioned

And also the below Matrix Identity

With

Take a note of the final result of the Gradient, which is the same form that we arrived at earlier under the Alternative-1 calculation of Gradient

And above is the exact formulae that we will implement in Python/Numpy very soon below.

So now let's go back to the original Cost Function

Which in Vectorized Form for the Mean Squared Error is defined as below

And after calculating the Gradient of this MSE in Vectorized form, which we did above the Gradient-Descent Algorithm will update the weights (θ / Theta values ) as below

Compare the above with the Gradient-Descent formulae for the Numerical case

A manual example of the Gradient-Descent implementation

Let's say for simple single variable training dataset we have the following values

x,y
1,1
2,2
3,3
4,4

Further, assume,

α (learning rate) = 1

m (number of training examples) =4

Setting θ0 to 0 and θ1 to 1

So we have the linear equation

Now by convention,

at the i−th row

Further I denote,

(i.e. after k repetitions of the GD algorithm).

So after a single update, with GD algorithm, i.e. applying the below

So, regardless of how many times I apply the GD algorithm, the value of θ1 will be constantly equal to 1, since at every iteration we have θ0=0 and θ1=1

Now extend the above to a multivariable case,

Assume theta values have been picked at random as below

And the training-dataset is

So here, first, to calculate the hypothesis Equation, I need to transpose θ to give our initial vector θ

And the GD-Algorithm is,

And for applying the GD algorithm again, I need to evaluate

Python/Numpy Implementation

Please refer to the jupyter notebook

First, generate a training dataset in Matrix form

NumPy zeros() function in above — you can create an array that only contains only zeros using the NumPy zeros() function with a specific shape. The shape is row by column format. Its syntax is as below

For example, the code to generate a Matrix of 2 by 3 (2 rows and 3 columns)

import numpy as npA = np.zeros(shape = (2, 3))
print# Output below
[[0. 0. 0.]
 [0. 0. 0.]]

Which produces an array like the following:

If I run the above gen_data() function above for a set of 5 training data-set as below with bias and variance of 20 and 10 respectively

gen_data(5, 20, 10)

I will have the following form of output

(array([[1., 0.],
 [1., 1.],
 [1., 2.],
 [1., 3.],
 [1., 4.]]), 
 array([22.38023816, 24.18406356, 28.01360908, 26.80051617, 29.30101971])
 )

And now the function for Gradient-Descent implementing the Grdient formulae for a Mactrix that we derived above

And now finally invoke the above 2 functions to create some linear data and run the gradient-descent and also plot it to a graph.

And the resultant graph of the Linear line will be this

The full notebook is here

Reference and Further Reading

Matrix Multiplication — https://en.wikipedia.org/wiki/Matrix_multiplication

Matrix-Calculus — https://en.wikipedia.org/wiki/Matrix_calculus

Vector_Field — https://en.wikipedia.org/wiki/Vector_field

Matrix Transpose Properties -https://en.wikipedia.org/wiki/Transpose#Properties

Matrix Cookbook — https://www.math.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf

Online Calculation of Matrix Derivative — http://www.matrixcalculus.org/