Chapter 3

Basic Mathematical Computations

This chapter discusses typical scientific computations using Python code. We will focus on how numerical data are represented mathematically, how they are structured or organized, stored, manipulated, processed, and computed in efficient ways. The subtleties of these operations in Python will be considered. All code presented in this textbook, and in this chapter in particular, can be found at https://sohoware.ru/SohoBook/. Our discussion begins with several basic linear-algebra operations on data structured as...

99min 18,652Words 15Materials

Key ideas

Linear Algebra
Scalar Numbers
Vectors
Matrices
Tensors
Dot Product of Two Vectors

Practice assignment

Take a small list or dictionary of data, convert it into a NumPy array, and compute two simple statistics. Connect the result with "Basic Mathematical Computations".

Open lab

Basic Mathematical Computations

Theme

Open lab

3.1 Linear Algebra

Linear algebra is extremely important for any computation involving large data, such as machine learning. We plan to briefly review the main linear-algebra operations using Python programming and modules that have already been developed by the broad Python community. We will go through the main concepts, mathematical notation, data structures, and computational procedures. Readers may skip this chapter if they are already confident in basic linear-algebra computations. Our discussion begins with data structures. First, we import the required modules and functions.

import numpy as np
# Import the NumPy package and assign it the alias np.

# We will also use the TensorFlow package below.
# If this has not yet been done, install TensorFlow with the command: pip install tensorflow.
# After installation, import TensorFlow.

import tensorflow as tf  # Import TensorFlow and give it the alias tf.

3.1.1 Scalar Numbers

In the previous chapter, we already touched on scalar numbers in the context of mathematical computations. These numbers are divided into three main categories: integers, real numbers, and complex numbers. In the Python programming language, each such number is assigned a unique name and memory address. This allows us to refer to the number by its name, update its value, and use it as an argument in functions, whether they are built-in functions, functions defined by the user in the code, or functions from imported modules.

Operations with numbers in Python are implemented in a way that is intuitive and straightforward. However, computations often encounter a problem related to exceeding the admissible numerical range, which can lead to overflow or underflow. This means that a number becomes too large or too small to be processed correctly under the limitations imposed by numerical types and the computer architecture.

In Python, numeric values are assumed to be able to cover the full range of real numbers, with precision limited by the characteristics of the computer. This is important in machine learning, where precision and the ability to process large numerical ranges may be essential for the efficiency and accuracy of models.

3.1.2 Vectors

A vector is understood as an ordered set of numbers arranged in one dimension. These numbers are called the components of the vector. In a physical context, each component of a vector may represent a quantity in a certain direction. For example, a force vector in three-dimensional space contains three components, each of which corresponds to the force acting along the respective coordinate axis.

In machine learning, vectors are used to represent feature sets, where each element of the vector corresponds to one feature. Such vectors can be very large, especially when methods such as the finite-element method are used, where objects are discretized into many elements and nodes, resulting in vectors with millions of components.

Let us consider vectors through a programming example using Python and its scientific-computing libraries, such as TensorFlow and NumPy.

Example of creating a vector in TensorFlow:

import tensorflow as tf  # Import TensorFlow and give it the alias tf.

p = 15  # Vector length, representing the number of components.
# Create a one-dimensional tensor from consecutive numbers from 0 to p - 1.
x = tf.range(p)  # Use the TensorFlow range function.
# Print the tensor and its shape.
print(x)  # Tensor x now contains 15 components.
print(x.shape)  # The tensor shape in Python is denoted as (p,).

tf.Tensor([ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14], shape=(15,), dtype=int32)
(15,)

As a result, we obtain a vector x of length 15. In mathematics, such a vector is often represented as a column, but in Python it is represented by default as a row.

Now, if we want to transpose the vector, that is, convert a row vector into a column vector, we first need to reshape it:

import tensorflow as tf  # Import TensorFlow and give it the alias tf.

# Change the tensor shape to make it two-dimensional, that is, a column vector.
x_reshaped = tf.reshape(x, (p, 1))
# Transpose the two-dimensional tensor.
x_transposed = tf.transpose(x_reshaped)
# Print the transposed tensor and its shape.
print(x_transposed)  # Transposed vector.
print(x_transposed.shape)  # Shape of the transposed vector.

tf.Tensor([[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14]], shape=(1, 15), dtype=int32)
(1, 15)

In TensorFlow, transposing a vector does not necessarily cause the data to be physically moved in memory, which makes the operation efficient. However, TensorFlow may not preserve identical memory addresses for the original and transposed tensors, especially when operations that change tensor shape are used.

3.1.3 Matrices

A matrix is a data structure consisting of numbers ordered along two dimensions: rows and columns. It is an extension of the vector concept to two or more dimensions. Each column of a matrix can be considered a vector. In this sense, a matrix can be viewed as a set of vectors stacked side by side. Matrices are used to represent and process datasets in which rows correspond to individual observations, for example different images, and columns correspond to different features, for example pixels.

In Python programming with the TensorFlow library, matrix operations can be performed as follows:

import tensorflow as tf  # Import TensorFlow and give it the alias tf.

p = 15  # Vector length, representing the number of components.
# Create a one-dimensional tensor from consecutive numbers from 0 to p - 1.
x = tf.range(p)  # Use the TensorFlow range function.
# Assuming that x is a one-dimensional tensor, reshape it
# to obtain a two-dimensional tensor with 3 rows and 5 columns.
A = tf.reshape(x, (3, 5))
# Print tensor A and its transposed version.
print('@' * 50)
print("A =", A)
print("A.T =", tf.transpose(A))
# Consider an example where we determine the shape of tensor A
# and of its transposed version.
print('@' * 50)
print(A.shape)
print(tf.transpose(A).shape)
# Tensor elements can be accessed
# using indexing and slicing:
print('@' * 50)
# Indexing the element at the intersection of the first row and the second column.
print('A[0, 1] =', A[0, 1])
print('row 2', A[1, :])  # Get all elements of the second row.
print('column 1', A[:, 1])  # Get all elements of the second column.

@@@@@@@@@@@@@@@@@@@@
A = tf.Tensor(
[[ 0  1  2  3  4]
 [ 5  6  7  8  9]
 [10 11 12 13 14]], shape=(3, 5), dtype=int32)
A.T = tf.Tensor(
[[ 0  5 10]
 [ 1  6 11]
 [ 2  7 12]
 [ 3  8 13]
 [ 4  9 14]], shape=(5, 3), dtype=int32)
@@@@@@@@@@@@@@@@@@@@
(3, 5)
(5, 3)
@@@@@@@@@@@@@@@@@@@@
A[0, 1] = tf.Tensor(1, shape=(), dtype=int32)
row 2 tf.Tensor([5 6 7 8 9], shape=(5,), dtype=int32)
column 1 tf.Tensor([ 1  6 11], shape=(3,), dtype=int32)

3.1.4 Tensors

In mathematics and physics, a tensor is a structured set of data that transforms according to certain rules when the coordinate system changes. Tensors may have different orders: a scalar, a single number, is a tensor of order zero; a vector is a tensor of the first order; a matrix is a tensor of the second order; and so on as the number of dimensions increases.

However, in the context of machine learning, a “tensor” usually denotes any data array whose dimensionality exceeds two. In this field, tensors are often associated with large data that must be structured in high dimension. For example, an RGB image is a three-dimensional tensor whose axes correspond to height, width, and color channels. In the Python NumPy library, tensors are represented as multidimensional arrays.

Machine learning usually does not use the tensor transformations that are applied in physics, so for the purposes of this textbook we will use the term “tensor” to denote multidimensional arrays without taking into account the mathematical rules of tensor transformations.

To create a tensor in the TensorFlow library, we can use the following code:

import tensorflow as tf  # Import TensorFlow and give it the alias tf.

# Create a three-dimensional tensor using tf.range, then reshape it.
X = tf.range(24)
X = tf.reshape(X, (2, 3, 4))
# Print the tensor dimension and the tensor itself.
print('X.shape =', X.shape)
print('X =', X)

X.shape = (2, 3, 4)
X = tf.Tensor(
[[[ 0  1  2  3]
  [ 4  5  6  7]
  [ 8  9 10 11]]

 [[12 13 14 15]
  [16 17 18 19]
  [20 21 22 23]]], shape=(2, 3, 4), dtype=int32)

The result is a three-dimensional tensor X with shape (2, 3, 4), representing two blocks, each of which contains three rows and four columns.

In machine learning, tensors are often used for operations applied element by element, as shown in the following example:

import tensorflow as tf  # Import TensorFlow and give it the alias tf.

# Create a two-dimensional tensor A and a tensor B
# containing ones of the same shape multiplied by 8.
A = tf.reshape(tf.range(8), (2, 4))
B = tf.ones_like(A) * 8
# Perform elementwise addition and multiplication of tensors A and B.
print('A =', A, '\nB =', B)
print('A + B =', A + B, '\nA * B =', A * B)

A = tf.Tensor(
[[0 1 2 3]
 [4 5 6 7]], shape=(2, 4), dtype=int32)
B = tf.Tensor(
[[8 8 8 8]
 [8 8 8 8]], shape=(2, 4), dtype=int32)
A + B = tf.Tensor(
[[ 8  9 10 11]
 [12 13 14 15]], shape=(2, 4), dtype=int32)
A * B = tf.Tensor(
[[ 0  8 16 24]
 [32 40 48 56]], shape=(2, 4), dtype=int32)

Here tensor A consists of the numbers from 0 to 7, while tensor B consists of the number 8 repeated in every element. Addition and multiplication are performed element by element, resulting in a new tensor with the same shape.

In machine learning, it is often necessary to compute the sum or mean value of elements in a tensor. In this case, a tensor can be treated as a multidimensional data array. Operations on tensors, such as summation and mean calculation, make it possible to obtain generalized characteristics of the data.

Let us consider how this operation is performed in practice using the TensorFlow library in Python:

import tensorflow as tf  # Import TensorFlow and give it the alias tf.

# Initialize a one-dimensional tensor, or vector, using the range function.
x = tf.range(5)
print(x)  # Print the vector.
print(tf.reduce_sum(x))  # Compute the sum of vector elements.

tf.Tensor([0 1 2 3 4], shape=(5,), dtype=int32)
tf.Tensor(10, shape=(), dtype=int32)

This code creates a one-dimensional tensor, or vector, x with values from 0 to 4. The output of print(x) shows this vector. The command print(tf.reduce_sum(x)) returns the sum of all elements of vector x, that is, 0 + 1 + 2 + 3 + 4, which is equal to 10.

Next, we create a two-dimensional tensor, a matrix, and perform operations:

import tensorflow as tf  # Import TensorFlow and give it the alias tf.

# Create a two-dimensional tensor filled with ones.
X = tf.ones((3, 5))
print(X)  # Print the matrix.
print(tf.reduce_sum(X))  # Compute the sum of matrix elements.
# Compute the mean value of matrix elements in two ways.
print(tf.reduce_mean(X), tf.reduce_sum(X) / tf.size(X).numpy())

tf.Tensor(
[[1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]], shape=(3, 5), dtype=float32)
tf.Tensor(15.0, shape=(), dtype=float32)
tf.Tensor(1.0, shape=(), dtype=float32) tf.Tensor(1.0, shape=(), dtype=float32)
<NDArray 1 @cpu(0)>
[1.]
<NDArray 1 @cpu(0)>

This code forms a 3 × 5 matrix X in which each element is equal to 1. The output of print(X) shows this matrix. The command print(tf.reducesum(X)) computes the sum of all elements of matrix X, which in this case is equal to 15 because the matrix contains fifteen ones. Here tf.reducemean(X) computes the mean value of all elements of tensor X, which in this case is 1 because all tensor elements are equal to one. The same result can be obtained by dividing the sum of the elements by the number of elements in the tensor: tf.reduce_sum(X) / tf.size(X).numpy().

These operations form the basis of data analysis, allowing the dimensionality of information to be reduced to manageable and interpretable statistical values.

3.1.5 Dot Product of Two Vectors

The dot product is the most frequently used operation in machine learning. In this section we will discuss in detail how it is used for vectors that may have different data structures, which introduces certain subtleties.

Given two vectors a and b, their dot product in linear algebra is often written as a · b. In essence, it is simply the sum of the products of their corresponding elements, which results in a scalar. This means that the shapes of vectors a and b must be compatible: they must have the same length. Let us consider several examples.

import numpy as np  # Import the NumPy package and assign it the alias np.
import tensorflow as tf  # Import TensorFlow and give it the alias tf.

# Initialize two vectors, a and b.
a = tf.range(5)
b = tf.ones_like(a) * 2  # This ensures their compatibility.
# Print the vectors and their shapes.
print('--------------------------------------')
print(f"a={a}, a.shape={a.shape}\nb={b}, b.shape={b.shape}")
# Compute the dot product of the vectors using the TensorFlow library.
print('--------------------------------------')
print('tf.tensordot(a, b, axes=1)', tf.tensordot(a, b, axes=1))
# Print the shape of the dot-product result.
print('--------------------------------------')
print('tf.tensordot(a, b, axes=1).shape', tf.tensordot(a, b, axes=1).shape)
# Compute the dot product of the vectors using the NumPy library.
print('--------------------------------------')
print(f"np.dot(a, b)={np.dot(a.numpy(), b.numpy())}")

--------------------------------------
a=[0 1 2 3 4], a.shape=(5,)
b=[2 2 2 2 2], b.shape=(5,)
--------------------------------------
tf.tensordot(a, b, axes=1) tf.Tensor(20, shape=(), dtype=int32)
--------------------------------------
tf.tensordot(a, b, axes=1).shape ()
--------------------------------------
np.dot(a, b)=20

The dot product of vectors a and b gives the same scalar result in TensorFlow and NumPy.

Transposition is a mathematical operation in which the rows of a matrix become columns, and vice versa. This operation is often used in mathematics and computer computations, especially when working with matrices.

However, when we are dealing with one-dimensional arrays, or vectors, in NumPy, the situation is slightly different. In NumPy, a one-dimensional array or vector has shape (n,), where n is the number of elements in the array. Only one dimension is specified, corresponding to the length of the vector. Since there is no second dimension, the transposition operation has no meaning because there is no second axis that could be “swapped” with the first. Transposition assumes a change in two-dimensional structure, but a one-dimensional array is already in its simplest form.

In a broader sense, when we transpose a two-dimensional array in NumPy, we swap its axes. If we have an array of shape (m, n), where m is the number of rows and n is the number of columns, transposition changes its shape to (n, m). This means that elements that were in row i and column j before transposition will be in row j and column i after transposition.

For a one-dimensional array, such a change of axes cannot be performed because there is only one dimension. The .T transposition property or the transpose() method in NumPy returns the original one-dimensional array unchanged because there is no axis to permute.

Here is a simple example for illustration:

import numpy as np

# Create a one-dimensional array in NumPy.
a = np.array([1, 2, 3])
# Try to transpose the one-dimensional array.
a_transposed = a.T
# The result of transposition is the same one-dimensional array.
print(a_transposed)  # Prints [1 2 3].

[1 2 3]

As you can see, the result of transposing a one-dimensional array is the same one-dimensional array, because it has only one axis and therefore no other axis with which it could be swapped. A one-dimensional NumPy array has shape (n,) and is not treated as a matrix.

When b is a column vector, a special case of a two-dimensional array, it has two axes, just like a matrix. The dot product a · b is the same as the matrix product ab, where a is defined as a row vector and b as a column vector, in terms of the resulting scalar value. Thus, in our presentation, we do not distinguish between them from a mathematical point of view and often use the equality a · b = ab.

Now let us provide code examples that illustrate the subtleties mentioned above:

import numpy as np  # Import the NumPy package and assign it the alias np.
import tensorflow as tf  # Import TensorFlow and give it the alias tf.

# Initialize two vectors, a and b.
a = tf.range(5)
b = tf.fill([5], 2)  # This ensures their compatibility.
# Convert the one-dimensional array b into a column vector.
b_c = tf.reshape(b, [-1, 1])
# Print the column vector and its shape.
print('-----------------------------------')
print(b_c, 'b_c.shape=', b_c.shape)
# Compute the dot product in NumPy; the result is a scalar.
print('-----------------------------------')
print('np.dot(a, b)', np.dot(a.numpy(), b.numpy()))
# Compute the dot product of a one-dimensional array and a column vector; the result is an array.
print('-----------------------------------')
print('np.dot(a, b_c)', np.dot(a.numpy(), b_c.numpy()))
# Print the results of dot-product operations in TensorFlow.
print('-----------------------------------')
print('tf.tensordot(a, b, axes=1)', tf.tensordot(a, b, axes=1))
print('-----------------------------------')
print('tf.tensordot(a, b_c, axes=1)', tf.tensordot(a, b_c, axes=1))
print('-----------------------------------')
print('tf.tensordot(a, b_c, axes=1).shape', tf.tensordot(a, b_c, axes=1).shape)

-----------------------------------
tf.Tensor(
[[2]
 [2]
 [2]
 [2]
 [2]], shape=(5, 1), dtype=int32) b_c.shape= (5, 1)
-----------------------------------
np.dot(a, b) 20
-----------------------------------
np.dot(a, b_c) [20]
-----------------------------------
tf.tensordot(a, b, axes=1) tf.Tensor(20, shape=(), dtype=int32)
-----------------------------------
tf.tensordot(a, b_c, axes=1) tf.Tensor([20], shape=(1,), dtype=int32)
-----------------------------------
tf.tensordot(a, b_c, axes=1).shape (1,)

As can be seen, all these operations give the same scalar value.

The dot product is a standard operation in linear algebra and has special significance when we work with vectors in column form. Vectors are usually represented as columns when we deal with linear algebra, and in this context the dot product of two column vectors results in a scalar. However, in mathematical notation this product is written using the transposition of one of the vectors to show that we are summing the products of corresponding elements of the two vectors.

Suppose we have two vectors, a and b, each of which has column shape and contains the same number of elements. The dot product of these vectors is denoted as aᵀb or bᵀa, where the symbol T denotes the transposition operation. Transposing vector a turns it from a column vector into a row vector, after which the products of corresponding elements of vector a and vector b are summed. This product gives us a scalar, that is, a single number.

It is important to note that although the result is a scalar value, in NumPy the result may be represented as a two-dimensional array or matrix with one element, which technically makes it a 1 × 1 matrix rather than just a number. This is because in these libraries vectors and matrices are represented as arrays, and operations on them follow array rules even when the result is a single value.

import numpy as np  # Import the NumPy package and assign it the alias np.
import tensorflow as tf  # Import TensorFlow and give it the alias tf.

# Initialize two vectors, a and b.
a = tf.range(5)
b = tf.fill([5], 2)  # This ensures their compatibility.
# Convert the one-dimensional array b into a column vector.
b_c = tf.reshape(b, [-1, 1])
# Convert the one-dimensional array a into a column vector.
a_c = tf.reshape(a, [-1, 1])
# Compute the dot product of the transposed column vector a
# and the column vector b; the result is a scalar in a 2D array.
print('-------------------------------------')
print('tf.matmul(a_c, b_c, transpose_a=True)', tf.matmul(a_c, b_c, transpose_a=True))
# To obtain a scalar value from the 2D array of shape (1, 1), use:
print('-------------------------------------')
print('scalar value', tf.reshape(tf.matmul(a_c, b_c, transpose_a=True), []).numpy())
# Use flattening to convert the column vector back into a row vector.
print('-------------------------------------')
print('tf.tensordot(a_c, b_c, axes=1)', tf.tensordot(tf.reshape(a_c, [-1]), tf.reshape(b_c, [-1]), axes=1))
# To obtain a scalar value from the one-dimensional array of shape (1,), use:
print('-------------------------------------')
print('scalar value', tf.tensordot(tf.reshape(a_c, [-1]), tf.reshape(b_c, [-1]), axes=1).numpy())

-------------------------------------
tf.matmul(a_c, b_c, transpose_a=True) tf.Tensor([[20]], shape=(1, 1), dtype=int32)
-------------------------------------
scalar value 20
-------------------------------------
tf.tensordot(a_c, b_c, axes=1) tf.Tensor(20, shape=(), dtype=int32)
-------------------------------------
scalar value 20

Thus, the code shown above demonstrates how the dot-product operation is performed in TensorFlow and NumPy, and how to handle the results of these operations in different data structures.

3.1.6 Outer Product of Two Vectors: Definition and Application

The outer product of two vectors a and b, denoted as a ⊗ b, results in the creation of a matrix. This operation is always compatible regardless of the shapes of vectors a and b.

How it works: the matrix element at position ij is the product ai bj. For example, if we have a vector with elements [0, 1, 2] and a vector with elements [2, 2, 2, 2, 2], the result of the outer product will be a matrix in which each element of one vector is multiplied by each element of the other.

Consider an example in Python using the NumPy library:

import numpy as np  # Import the NumPy package and assign it the alias np.

# Define vectors a and b.
a = np.arange(3)        # Creates vector [0, 1, 2].
b = np.ones(5) * 2      # Creates vector [2, 2, 2, 2, 2].
# Print vectors a and b.
print('-------------------------------')
print(a, b)
# Compute the outer product of a and b.
# The np.outer() function takes two vectors and returns their outer product.
outer_product = np.outer(a, b)
print('-------------------------------')
print('np.outer=\n', outer_product)

-------------------------------
[0 1 2] [2. 2. 2. 2. 2.]
-------------------------------
np.outer=
 [[0. 0. 0. 0. 0.]
 [2. 2. 2. 2. 2.]
 [4. 4. 4. 4. 4.]]

A similar result can also be achieved using the @ operator by converting vector a into a column, shape (n, 1), and vector b into a row, shape (1, m).

import numpy as np  # Import the NumPy package and assign it the alias np.

# Define vectors a and b.
a = np.arange(3)        # Creates vector [0, 1, 2].
b = np.ones(5) * 2      # Creates vector [2, 2, 2, 2, 2].
# Convert the one-dimensional array a into a column vector.
a_c = a.reshape(-1, 1)
# Convert the one-dimensional array b into a column vector.
b_c = b.reshape(-1, 1)
# Transpose the column vector b.
b_transposed = b_c.T
# Print vectors a_c and b_transposed.
print('-------------------------------')
print(a_c, b_transposed)
# Compute the outer product of a and b.
# result is the outer-product result.
result = a_c @ b_transposed
print('-------------------------------')
print('np.outer=\n', result)

-------------------------------
[[0]
 [1]
 [2]] [[2. 2. 2. 2. 2.]]
-------------------------------
np.outer=
 [[0. 0. 0. 0. 0.]
 [2. 2. 2. 2. 2.]
 [4. 4. 4. 4. 4.]]

However, using the built-in np.outer() function is preferable because it works faster and does not require additional operations.

3.1.7 Matrix-Vector Product

The product of a matrix and a vector is a standard operation in linear algebra that is performed when the dimensions of the matrix and the vector are compatible. This operation makes it possible to transform data by multiplying them by feature weights represented as a vector.

Let us consider how this operation is implemented in practice using the TensorFlow library in Python:

import tensorflow as tf  # Import TensorFlow and give it the alias tf.

# Initialize matrix A35 of size 3x5 with the numbers from 0 to 14.
A35 = tf.reshape(tf.range(15), (3, 5))
# Initialize vector d5 of length 5 consisting of ones.
d5 = tf.cast(tf.ones((A35.shape[1],)), tf.int32)  # Convert d5 to int32.
# Print matrix A35, its shape, vector d5, and its shape.
print('--------------------')
print(A35.numpy(), A35.shape, d5.numpy(), d5.shape)
# Compute the product of matrix A35 and vector d5.
f = tf.matmul(A35, tf.reshape(d5, [-1, 1]))
# Shapes are compatible: [3, 5] x [5, 1] -> vector of length 3.
print('--------------------')
print(f.numpy(), f.shape)

--------------------
[[ 0  1  2  3  4]
 [ 5  6  7  8  9]
 [10 11 12 13 14]] (3, 5) [1 1 1 1 1] (5,)
--------------------
[[10]
 [35]
 [60]] (3, 1)

In this example, we multiply matrix A35 of size 3 × 5 by vector d5 of length 5. Since the number of columns in matrix A35 corresponds to the length of vector d5, the multiplication operation can be performed, and the result is a vector of length 3.

However, if we try to perform the operation in the other order, that is, multiply vector d5 by matrix A35, we obtain an error because the dimensions are not compatible for such an operation.

Replace the variable f:

# The following operation will produce an error because of incompatible shapes.
# f = tf.matmul(d5, A35)  # This line will raise an error.

If we want to multiply a vector on the left by a matrix, the vector must have a length equal to the number of rows of the matrix:

import tensorflow as tf  # Import TensorFlow and give it the alias tf.

# Initialize matrix A35 of size 3x5 with the numbers from 0 to 14.
A35 = tf.reshape(tf.range(15), (3, 5))
# Initialize vector d3 of length 3 consisting of ones.
# Use the first element of the A35 shape to determine the vector length.
d3 = tf.cast(tf.ones((A35.shape[0],)), tf.int32)
# Print matrix A35, its shape, vector d3, and its shape.
print('--------------------')
print(d3, d3.shape, A35, A35.shape)
# Compute the product of vector d3 and matrix A35.
f = tf.matmul(tf.reshape(d3, [1, -1]), A35)
# Shapes are compatible: [1, 3] x [3, 5] -> vector of length 5.
print('--------------------')
print(f, f.shape)

--------------------
tf.Tensor([1 1 1], shape=(3,), dtype=int32) (3,) tf.Tensor(
[[ 0  1  2  3  4]
 [ 5  6  7  8  9]
 [10 11 12 13 14]], shape=(3, 5), dtype=int32) (3, 5)
--------------------
tf.Tensor([[15 18 21 24 27]], shape=(1, 5), dtype=int32) (1, 5)

In this case, the result is a vector of length 5, obtained by summing the elements of each column of matrix A35.

3.1.8 Matrix-Matrix Multiplication

To multiply two matrices, their shapes must be compatible. This means that the number of columns in the first matrix must match the number of rows in the second matrix.

Consider an example of multiplying a matrix by a matrix using the TensorFlow library in Python:

import tensorflow as tf  # Import TensorFlow and give it the alias tf.

# Initialize two matrices, A23 and B35, with compatible shapes.
A23 = tf.ones(shape=(2, 3))
B35 = tf.ones(shape=(3, 5))
print('-' * 20)
# Print both matrices.
print(A23.numpy(), B35.numpy())
# Perform the matrix-by-matrix multiplication operation.
C25 = tf.matmul(A23, B35)
# Shapes are compatible: [2, 3] x [3, 5] -> result has shape [2, 5].
print('-' * 20)
print(C25.numpy())

--------------------
[[1. 1. 1.]
 [1. 1. 1.]] [[1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]]
--------------------
[[3. 3. 3. 3. 3.]
 [3. 3. 3. 3. 3.]]

In this code, matrix A23 has shape 2 × 3, and matrix B35 has shape 3 × 5. After multiplication we obtain matrix C25 of shape 2 × 5, where each element is the sum of products of elements of the corresponding rows of the first matrix and columns of the second matrix. In this case, because both matrices are filled with ones, each element of the resulting matrix C25 is equal to 3.

The same operation in NumPy is performed similarly:

import numpy as np  # Import the NumPy package and assign it the alias np.

# Initialize two matrices, A23 and B35, with compatible shapes.
A23 = np.ones(shape=(2, 3))
B35 = np.ones(shape=(3, 5))
print('-' * 20)
# Print both matrices.
print(A23, '\n', B35)

print('-' * 20)
# Use the np.dot() function to multiply matrices.
print('np.dot():\n', np.dot(A23, B35))
print('-' * 20)
# Use the @ operator in NumPy to multiply matrices.
print('numpy @ operator:\n', A23 @ B35)

--------------------
[[1. 1. 1.]
 [1. 1. 1.]]
 [[1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]]
--------------------
np.dot():
 [[3. 3. 3. 3. 3.]
 [3. 3. 3. 3. 3.]]
--------------------
numpy @ operator:
 [[3. 3. 3. 3. 3.]
 [3. 3. 3. 3. 3.]]

Both approaches in NumPy produce the same result as in TensorFlow.

3.1.9 Norms

In linear algebra and data analysis, a norm is a function that assigns a positive length or magnitude to every vector in a vector space, except for the zero vector, which is assigned length zero. In the context of Python programming and machine learning, norms are used to determine the “magnitude” of vectors and matrices. They help determine how large the elements of a vector or matrix are relative to zero or to one another.

There are several different types of norms, each of which has its own applications.

The L2 norm, also known as the Euclidean norm, of a vector is the square root of the sum of squares of its elements. It is the most frequently used type of norm, especially in machine learning, because it is related to the Euclidean distance between points. For matrices, the L2 norm is also known as the Frobenius norm, which is the analogue of the L2 norm for matrix data.

The L1 norm of a vector is the sum of the absolute values of its elements. This norm is less sensitive to outliers and is often used in statistics and optimization. For matrices, the L1 norm can be defined as the maximum L1 norm among the columns of the matrix.

Example of using the L1 norm and L2 norm in Python with TensorFlow:

import tensorflow as tf  # Import TensorFlow and give it the alias tf.
import numpy as np  # Import the NumPy package and assign it the alias np.

# Create a tensor, or vector, of ones.
d = tf.ones(9)
print('-' * 20)
# Print the vector and the sum of its elements.
print(d.numpy(), tf.reduce_sum(d).numpy())
print('-' * 20)
# Compute the L2 norm of the vector using TensorFlow.
print('tf.norm(d)', tf.norm(d).numpy())
print('-' * 20)
# Compute the L2 norm of the vector using NumPy.
print('np.linalg.norm(d)', np.linalg.norm(d.numpy()))
print('-' * 20)
# Compute the L1 norm of the vector using TensorFlow.
print('tf.norm(d, ord=1)', tf.norm(d, ord=1).numpy())
print('-' * 20)
# Compute the L1 norm of the vector using NumPy.
print('np.linalg.norm(d, 1)', np.linalg.norm(d.numpy(), 1))

--------------------
[1. 1. 1. 1. 1. 1. 1. 1. 1.] 9.0
--------------------
tf.norm(d) 3.0
--------------------
np.linalg.norm(d) 3.0
--------------------
tf.norm(d, ord=1) 9.0
--------------------
np.linalg.norm(d, 1) 9.0

Norms are used in regularization, which prevents overfitting of machine-learning models by limiting the magnitude of model parameters. For example, L1 regularization can lead to sparse solutions, where some model parameters become equal to zero; this may be useful for feature selection. L2 regularization tends to reduce all model parameters simultaneously, which helps avoid excessively large weights and reduces model complexity.

3.1.10 Solving Systems of Linear Algebraic Equations

In mathematics and machine learning, the problem of solving systems of linear equations often arises. Such systems can be represented in the form of the matrix equation:

KD = F

Here K is a square matrix, D is the unknown vector of variables, and F is the vector of results or known values.

Matrix K must be positive definite, which means that all its eigenvalues are greater than zero. This condition is required to ensure uniqueness of the solution. In the context of the finite-element method, or FEM, the stiffness matrix K is usually symmetric and positive definite, which guarantees that the equations of the system correspond to a physically real and stable model.

The function numpy.linalg.solve is intended for solving such systems of equations. It takes two arguments, matrix K and vector F, and returns vector D, which satisfies the original equation.

Let us consider an example. Suppose we have a system of equations where matrix K has certain physical characteristics of the system, and F represents external forces acting on the system. We want to find D.

import numpy as np  # Import the NumPy package and assign it the alias np.

# Define a new matrix K.
K = np.array([[3, 1], [1, 2]])
# Define vector F.
F = np.array([9, 3])
# Solve the system of equations and find D.
D = np.linalg.solve(K, F)
# Print the results to the screen.
print("Matrix K:\n", K)
print("Vector F:", F)
print("Vector D:", D)

Stiffness matrix K:
 [[3 1]
 [1 2]]
External-load vector F: [9 3]
Vector D: [3. 0.]

In this code, numpy.linalg.solve efficiently solves the system of equations. If matrix K were not square, or if the system were overdetermined, that is, had more equations than unknowns, we would use numpy.linalg.lstsq to find a least-squares solution that minimizes the sum of squared errors between data and model predictions. This approach is often used in machine learning for regression problems.

It is worth noting that solving large systems of equations can be a resource-intensive process. Modern numerical algorithms and advances in computing technology make it possible to handle this task efficiently, especially when iterative methods and gradient-descent-based methods are used to reduce the solution error to an acceptable level. These methods are similar to those used in machine learning to optimize loss functions.

3.1.11 Matrix Inversion

Matrix inversion is the process of finding a matrix K⁻¹ that is inverse to a given matrix K. The inverse matrix is unique for K and satisfies the conditions K × K⁻¹ = I and K⁻¹ × K = I, where I is the identity matrix of the corresponding size.

In the previous sections we discussed the solution of linear-algebra equations. Let us consider matrix inversion in more detail. Consider the equation KD = F, whose solution can be represented as:

D = K⁻¹F,

where K⁻¹ is the inverse matrix of K. Thus, if we can compute K⁻¹, the solution is reduced to the product of a matrix and a vector. For small systems this approach indeed works and is widely used. To compute the inverse matrix we use the numpy.linalg.inv() function:

import numpy as np  # Import the NumPy package and assign it the alias np.
from numpy.linalg import inv

# Define matrix K.
K = np.array([[3, 1], [1, 2]])
# Compute the inverse matrix of K.
Kinv = inv(K)
print('-' * 20)
print(Kinv)
print('-' * 20)
# Check that the product of K and Kinv gives the identity matrix.
print(np.allclose(np.dot(K, Kinv), np.eye(2)))
print(np.allclose(np.dot(Kinv, K), np.eye(2)))

--------------------
[[ 0.4 -0.2]
 [-0.2  0.6]]
--------------------
True
True

Here np.allclose is used to check that the product of K and Kinv gives approximately the identity matrix, which is a sign of correct inversion.

The solution of the equation can be obtained as follows:

import numpy as np  # Import the NumPy package and assign it the alias np.
from numpy.linalg import inv

# Define matrix K.
K = np.array([[3, 1], [1, 2]])
# Define vector F.
F = np.array([9, 3])
# Compute the inverse matrix of K.
Kinv = inv(K)
# Print result D to the screen.
D = np.dot(Kinv, F)
print('D:', D)

D: [3. 0.]

This solution matches the one obtained earlier. An interesting point is that several matrices can be inverted at the same time:

import numpy as np  # Import the NumPy package and assign it the alias np.
from numpy.linalg import inv

a = np.array([[[1., 2.], [3., 4.]], [[5, 6], [7, 8]]])
print('-' * 20)
print(a)
print('-' * 20)
print(inv(a))

--------------------
[[[1. 2.]
  [3. 4.]]

 [[5. 6.]
  [7. 8.]]]
--------------------
[[[-2.   1. ]
  [ 1.5 -0.5]]

 [[-4.   3. ]
  [ 3.5 -2.5]]]

Inverting large matrices is very expensive in terms of computational resources, so it is often preferable to solve algebraic systems of equations directly. In computations related to machine learning, one may also encounter singular matrices that have no inverse, which causes computation to fail. Matrices are often nearly singular; this may allow computations to continue, but may also lead to serious errors that appear as unexpected or strange behavior. When such behavior is observed, there is a high probability that the system matrix may be “ill-conditioned,” and one should check for possible errors in the data or formulation procedure that may have led to a nearly singular system matrix. If the problem lies in the data themselves, data cleaning or checking for errors in the data may be required. After that, mathematical methods can be used, one of which is singular value decomposition, or SVD, to obtain the best possible information from the data.

The main point we want to emphasize here is that the most important factor determining whether a system of algebraic equations can be solved and whether a high-quality solution can be obtained is the property, characteristic, or condition of the system matrix. Eigenvalues, if they exist, and the corresponding eigenvectors are characteristics of the matrix.

3.1.12 Eigenvalue Decomposition of a Matrix

Basics of decomposition

Decomposition of a matrix into eigenvalues and eigenvectors is a process in which a diagonalizable matrix is represented as the product of a matrix of eigenvectors, a diagonal matrix of eigenvalues, and the transposed matrix of eigenvectors. This method is especially effective for real symmetric matrices, because in this case the eigenvalues are real numbers and the eigenvectors can be orthonormal.

Matrix of eigenvectors (V). The eigenvectors of a matrix represent directions in which the action of that matrix can be expressed as simple scaling. If a matrix A has an eigenvector v, then when A is multiplied by v, the result is v multiplied by some scalar value. This scalar value is known as the eigenvalue. V consists of the eigenvectors of matrix A arranged as its columns. For a symmetric matrix these vectors will be orthonormal, meaning that they have unit length and are orthogonal to each other.

Diagonal matrix of eigenvalues (Λ). This component of the decomposition is a diagonal matrix in which each element on the main diagonal is an eigenvalue corresponding to an eigenvector from matrix V. Eigenvalues reflect the factor by which an eigenvector is scaled when multiplied by matrix A. In the case of a symmetric matrix, all eigenvalues are real numbers.

Transposed matrix of eigenvectors (Vᵀ). The last component of the decomposition is the transposed matrix V. Transposing a matrix means replacing its rows with columns, or vice versa. In the case of an orthonormal matrix V for a symmetric matrix A, transposition is equivalent to finding the inverse matrix of V. This is because multiplying V by Vᵀ, or vice versa, gives the identity matrix.

For a symmetric positive-definite matrix A, the decomposition has the following form:

A = VΛVᵀ,

where V is the matrix of eigenvectors, Λ is the diagonal matrix of eigenvalues, and Vᵀ is the transposed matrix V.

Since V is orthonormal, the following holds:

VᵀV = I,

and, consequently, the inverse matrix of V is equal to its transpose:

V⁻¹ = Vᵀ.

Computing the inverse matrix

After matrix decomposition, computing the inverse becomes simple. From the definition of the inverse matrix, AA⁻¹ = I, and using the equations A = VΛVᵀ and V⁻¹ = Vᵀ, we obtain:

A⁻¹ = VΛ⁻¹Vᵀ.

The inverse matrix of the diagonal matrix Λ consists of the reciprocals of its diagonal elements.

Let us consider examples of computing eigenvalues and eigenvectors in Python:

import numpy as np  # Import the NumPy library.
from numpy import linalg as lg  # Import the linalg module from NumPy under the alias lg.

A = np.array([[1, 0.4, 0.8], [0.4, 1, 0.5], [0.8, 0.5, 1]])  # Create a 3x3 matrix.
# Call the eig function from the linalg module, which computes the eigenvalues (e)
# and eigenvectors (v) of matrix A.
e, v = lg.eig(A)
print('-' * 20)
# Print matrix A.
print('matrix A:\n', A)
print('-' * 20)
# Print the eigenvalues of matrix A.
print('Eigenvalues:', e)
print('-' * 20)
# Print the eigenvectors of matrix A.
print('Eigenvectors:\n', v)

--------------------
matrix A:
 [[1.  0.4 0.8]
 [0.4 1.  0.5]
 [0.8 0.5 1. ]]
--------------------
Eigenvalues: [2.15226471 0.19102751 0.65670778]
--------------------
Eigenvectors:
 [[-0.60615713 -0.66545116  0.43560106]
 [-0.48421405 -0.12572504 -0.86586949]
 [-0.63095982  0.73577712  0.24601167]]

Here we obtain three real eigenvalues and orthonormal eigenvectors.

Using the eigenvalues and eigenvectors, the original matrix can be reconstructed:

import numpy as np  # Import the NumPy library.
from numpy import linalg as lg  # Import the linalg module from NumPy under the alias lg.

A = np.array([[1, 0.4, 0.8], [0.4, 1, 0.5], [0.8, 0.5, 1]])  # Create a 3x3 matrix.
# Call the eig function from the linalg module, which computes the eigenvalues (e)
# and eigenvectors (v) of matrix A.
e, v = lg.eig(A)
# Create diagonal matrix lamd from eigenvalues e. The elements of lamd are on the main diagonal,
# while all other elements are zero.
lamd = np.diag(e)
# Reconstruct the original matrix A using the decomposition into eigenvalues and eigenvectors.
# Multiply the matrix of eigenvectors v by the diagonal matrix of eigenvalues lamd and by the transposed
# matrix of eigenvectors v.T. The @ operator denotes matrix multiplication.
A_recovered = v @ lamd @ v.T
# Print the reconstructed matrix A.
print('Reconstructed matrix A:\n', A_recovered)

Reconstructed matrix A:
 [[1.  0.4 0.8]
 [0.4 1.  0.5]
 [0.8 0.5 1. ]]

This demonstrates that the matrix reconstructed in this way coincides with the result obtained using the original matrix representation. The same method can also be applied to asymmetric matrices, but in that case the eigenvalues may be complex.

Decomposition into eigenvalues and eigenvectors is a powerful tool, especially when it comes to matrix analysis and transformation. It is one form of matrix decomposition and has particular significance for symmetric and positive-definite matrices.

3.1.13 Matrix Condition Number

The condition number is a measure that estimates how “sensitive” the results of a system of equations are to small changes in the input data. In the context of a matrix, this number shows how much the solution of a system of equations changes if small changes are made to the coefficients of the equations, that is, to the elements of the matrix. In real computations, rounding errors and inaccuracies in data are common. The condition number helps us understand how strongly these small errors will affect the final result. A small condition number means that even with errors the results will be fairly accurate. A high condition number indicates that the system of equations or mathematical problem is numerically unstable. This means that even small changes in the input data can lead to significant changes in the results.

Imagine a matrix as a set of linear equations. If you have a system of linear equations, you use a matrix to represent and solve them. When you solve this system, you expect to obtain certain values of the variables. If the matrix has a low condition number, it means that small changes in the coefficients, or matrix elements, will lead to small changes in the solution. For example, if you have a matrix with condition number 2, this means that if you change the input data by 1%, the results will change by roughly 2%. Conversely, if a matrix has a high condition number, small changes in the coefficients may lead to very large and unpredictable changes in the solution. For example, if a matrix has condition number 10,000, then a 1% change in the input data may change the results by 10,000%.

Examples in Python using NumPy:

import numpy as np  # Import the NumPy package and assign it the alias np.
from numpy import linalg as lg

A = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
print('-' * 20)
# Call the eig function from the linalg module, which computes the eigenvalues (e)
# and eigenvectors (v) of matrix A.
e, v = lg.eig(A)
print('Eigenvalues:', e)
print('-' * 20)
# Identity matrix.
# Obviously, it has three eigenvalues equal to 1.0.
print(A, '\n Condition number of A=', lg.cond(A))

--------------------
Eigenvalues: [1. 1. 1.]
--------------------
[[1 0 0]
 [0 1 0]
 [0 0 1]]
 Condition number of A= 1.0

Since matrix A is the identity matrix, we obtain condition number 1, as expected. Another example:

import numpy as np  # Import the NumPy package and assign it the alias np.
from numpy import linalg as lg

A = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 0]])
print('-' * 20)
# Call the eig function from the linalg module, which computes the eigenvalues (e)
# and eigenvectors (v) of matrix A.
e, v = lg.eig(A)
print('Eigenvalues:', e)
print('-' * 20)
print(A, '\nCondition number of A=', lg.cond(A))

--------------------
Eigenvalues: [1. 1. 0.]
--------------------
[[1 0 0]
 [0 1 0]
 [0 0 0]]
 Condition number of A= inf

Because matrix A is singular, its condition number is inf, which in NumPy means infinity, as expected.

import numpy as np  # Import the NumPy package and assign it the alias np.
from numpy import linalg as lg

A = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 10]])
print('-' * 20)
# Call the eig function from the linalg module, which computes the eigenvalues (e)
# and eigenvectors (v) of matrix A.
e, v = lg.eig(A)
print('Eigenvalues:', e)
print('-' * 20)
# Identity-like diagonal matrix with one larger entry.
print(A, '\n Condition number of A=', lg.cond(A))

--------------------
Eigenvalues: [ 1.  1. 10.]
--------------------
[[ 1  0  0]
 [ 0  1  0]
 [ 0  0 10]]
 Condition number of A= 10.0

The condition number of this matrix A is 10.0, corresponding to 10.0 / 1.0. Again, the condition number is the ratio of the largest eigenvalue to the smallest. The number 10 indicates that the matrix is less stable than the identity matrix but is not singular. We can conclude that if the largest eigenvalue of a matrix is very large or the smallest eigenvalue of a matrix is very small, the matrix is likely to be singular depending on their ratio.

This conclusion implies that normalizing a matrix, which is often used in machine learning, theoretically does not change its condition number. It may help reduce the loss of significant digits caused by the limited representation of floating-point numbers in computer hardware.

3.1.14 Matrix Rank

The rank of a matrix is a measure that determines the maximum number of linearly independent rows or columns in the matrix. Linearly independent rows or columns are those that cannot be expressed as a linear combination of other rows or columns. If a square matrix has full rank, this means that all its columns, or rows, are mutually linearly independent. Such a matrix is not singular. For a singular matrix, that is, when it has no inverse matrix, the rank must be less than the full rank. In the context of linear algebra, this means that the rows or columns of the matrix are not linearly independent. The difference between full rank, which is the maximum possible number of linearly independent rows or columns, and the actual rank is called rank deficiency. Imagine a 3 × 3 matrix that has only two linearly independent rows or columns. Such a matrix has rank 2. Since the full rank for a square 3 × 3 matrix must be 3 and the actual rank of the matrix is 2, the rank deficiency is 3 - 2 = 1.

For a non-square matrix, full rank is the number of its columns or rows, whichever is smaller. Similarly, it can also be rank-deficient if its rank is less than full rank.

Now let us examine this in more detail using NumPy.

import numpy as np  # Import the NumPy package and assign it the alias np.
from numpy.linalg import matrix_rank, eig

A = np.eye(4)  # Identity matrix.
rank = matrix_rank(A)  # Compute the rank.
print('Rank of the identity matrix=', rank)

Rank of the identity matrix= 4

Here we can see that the identity matrix has size 4 × 4 and full rank.

import numpy as np  # Import the NumPy package and assign it the alias np.
from numpy.linalg import matrix_rank, eig

A = np.array([[1, -0.2, 0], [0.1, 1, -0.5], [0.1, 1, -0.5]])  # Singular matrix.
rank = matrix_rank(A)  # Compute the rank.
print('Rank of the singular matrix=', rank)

Rank of the singular matrix= 2

This singular matrix has two linearly independent columns and therefore rank 2. It has rank deficiency 1. Consequently, it must also have one zero eigenvalue, as shown below. If a matrix has rank deficiency n, it will have n zero eigenvalues. This is easy to check using NumPy.

import numpy as np  # Import the NumPy package and assign it the alias np.
from numpy.linalg import matrix_rank
from numpy import linalg as lg

# Singular matrix.
A = np.array([[1, -0.2, 0], [0.1, 1, -0.5], [0.1, 1, -0.5]])
# Compute the rank.
rank = matrix_rank(A)
# Compute the eigenvalues.
e, v = lg.eig(A)
# Print values with 4 digits after the decimal point.
formatted_eigenvalues = np.around(e, 4)
print('eigenvalues of matrix A:', formatted_eigenvalues)
print('Rank of the singular matrix=', rank)

eigenvalues of matrix A: [ 0.9562  0.5438 -0.    ]
Rank of the singular matrix= 2

Example of a non-square matrix:

import numpy as np  # Import the NumPy package and assign it the alias np.
from numpy.linalg import matrix_rank

# Non-square matrix.
A = np.array([[1, -0.2, 2], [0.1, 1, -0.5]])
# Compute the rank.
rank = matrix_rank(A)
print('Rank of the non-square matrix=', rank)

Rank of the singular matrix= 2

This matrix has only two rows and rank 2. It has full rank.

3.1.15 Rotation Matrix

In the context of a two-dimensional coordinate system, a rotation matrix is used to change the position of points or vectors. This process is important, for example, when rotating images or graphs. The rotation matrix T is represented as:

T = [[cos θ, -sin θ],
     [sin θ,  cos θ]]

Here θ denotes the rotation angle. This angle determines how far and in which direction we want to rotate our object, for example a vector. The matrix consists of the cosine and sine of this angle arranged in such a way that it rotates a vector by the specified angle in two-dimensional space.

The vector d that we want to rotate is represented by two components in the coordinate system:

d = [u, v]

Now let us consider how this is implemented in Python using the NumPy library.

import numpy as np  # Import the NumPy package and assign it the alias np.

theta = 45  # Rotation angle in degrees. Here it is 45 degrees.
# Convert the angle from degrees to radians,
# because NumPy cos and sin functions work with radians.
theta_rad = np.deg2rad(theta)
# Compute the cosine and sine of the angle for the matrix elements.
cos, sin = np.cos(theta_rad), np.sin(theta_rad)
# Create a rotation matrix using cosine and sine.
T = np.array([[cos, -sin], [sin, cos]])
print('-' * 20)
print('Transformation matrix T:\n', T)
d = np.array([1, 0])  # Initial vector directed along the X axis.
print('-' * 20)
print('Initial vector d:', d)
print('-' * 20)
K = T @ d  # Rotate vector d by angle theta.
print('vector d rotated by angle theta\n', K)
print('-' * 20)
K = T @ (T @ d)  # Rotate vector d by angle 2 * theta.
print('vector d rotated by angle 2 * theta\n', K)
print('-' * 20)
T_2 = T @ T  # Combine two rotation matrices, which is
K = T_2 @ d  # equivalent to writing T @ (T @ d).
print('vector d rotated by angle 2 * theta, second method\n', K)

--------------------
Transformation matrix T:
[[ 0.70710678 -0.70710678]
 [ 0.70710678  0.70710678]]
--------------------
Initial vector d: [1 0]
--------------------
vector d rotated by angle theta
 [0.70710678 0.70710678]
--------------------
vector d rotated by angle 2 * theta
 [0. 1.]
--------------------
vector d rotated by angle 2 * theta, second method
 [0. 1.]

This code fragment demonstrates how sequential application of a rotation matrix leads to an increase in the rotation angle. In this case, two applications of the rotation matrix rotate the vector by an angle twice as large as the original.

An interesting point, however, is that if the angle θ is 45 degrees, then a rotation by 8 × 45 = 360 degrees, written as T @ (T @ (T @ (T @ (T @ (T @ (T @ T)))))), effectively returns the vector to its original position, because a rotation by 360 degrees does not change the position of the object.

These code examples show how a rotation matrix can be used to change the orientation of vectors in two-dimensional space. This tool is used in machine learning and computer vision, where such transformations can be used for data preparation and augmentation, as well as to improve the perception and analysis of images.

3.2 Interpolation

Interpolation is a commonly used numerical technique that makes it possible to obtain approximate values on the basis of known data. In a certain sense, machine learning is similar to interpolation. This section considers general issues related to interpolation using the NumPy library. Interpolation is also known as curve fitting. Here we show several examples of function interpolation and approximation using values specified at discrete points in space.

The first example uses numpy.interp.

import numpy as np  # import the NumPy package and assign it the alias np

# Data points for which the values are known
x_known = [1, 2, 3, 4, 5]
y_known = [2.1, 2.9, 3.8, 5.1, 7.2]

# Points for which we want to interpolate values
x_new = [1.5, 2.5, 3.5]

# Use numpy.interp for interpolation
y_new = np.interp(x_new, x_known, y_known)
print(y_new)

[2.5  3.35 4.45]

In this code, xknown and yknown are the known data. In this case they define a function whose values are known at discrete points. xnew contains new points for which we want to compute values by using interpolation. The function np.interp takes the new points xnew, together with the arrays xknown and yknown, and computes the interpolated values y_new.

This simple example demonstrates the basic concept of interpolation: finding values at new points on the basis of known values at other points. This is used in machine learning and data analysis when it is necessary to work with incomplete or discrete data.

3.2.1 Piecewise-Linear Interpolation

Piecewise-linear interpolation is a method in mathematics and numerical analysis that is used to find new points on the basis of a discrete set of known points. In the one-dimensional case, this means that we have several known points on a line and want to estimate values at intermediate points.

import numpy as np  # import the NumPy package and assign it the alias np

# Available data
xn = [1, 15, 50]  # data: given x coordinates
fn = [1, 20, 40]  # data: given function values at x points

# Query/predict f at a new point x
x = 25
f = np.interp(x, xn, fn)  # obtain the approximate value at point x
print(f'f({x:.3f})≈{f:.3f}')  # print the result
print('-' * 20)

x_2 = [11, 21, 31, 19, 10]
f_2 = np.interp(x_2, xn, fn)  # query at a larger number of points
print(f_2)

Textbook illustration: Basic Mathematical Computations — 3.2.1 Piecewise-Linear Interpolation

f(25.000)≈25.714
--------------------
[14.57142857 23.42857143 29.14285714 22.28571429 13.21428571]

In practice, interpolation can be a rather risky operation, so special care is required, especially during extrapolation.

Extrapolation is the process of extending function estimates beyond the range of known values. Unlike interpolation, extrapolation assumes the prediction of function values outside the given data points.

To avoid extrapolation, or at least be aware of it, a warning can be set when interpolation occurs outside the domain covered by the data.

import numpy as np  # import the NumPy package and assign it the alias np

# Available data
xn = [1, 15, 50]  # data: given x coordinates
fn = [1, 20, 40]  # data: given function values at x points

out_of_domain = -1111111.1  # a number used as a warning
print(np.interp(29, xn, fn, right=out_of_domain))  # print a number for extrapolation
print(np.interp(60, xn, fn, right=out_of_domain))  # warning

28.0
-1111111.1

Interpolation using higher-order polynomials may be more accurate, but it may also create more serious difficulties. Piecewise-linear approximation is often much safer and can be very effective when “dense data” are available. Below is an example of using piecewise-linear interpolation to approximate a sinusoidal function.

import numpy as np  # import the NumPy package and assign it the alias np
import matplotlib.pyplot as plt  # module for displaying results

x = np.linspace(0, 4 * np.pi, 40)  # data: x values
y = np.sin(x)  # data: function values at x points

# generate dense x data whose values are obtained through interpolation
xvals = np.linspace(0, 4 * np.pi, 50)
yinterp = np.interp(xvals, x, y)

plt.plot(x, y, 'o')  # display the original data points
plt.plot(xvals, yinterp, '-x')  # display interpolated data points
plt.show()  # show plots

Here x is an array of values from 0 to 4*pi, divided into 40 equal parts. y contains the sine-function values at these points. Then xvals is generated as a denser set of points, and np.interp is used to interpolate sine values at these new points. After interpolation, the results are visualized with matplotlib. The original data points are shown as circles ('o'), and the interpolated points are shown as crosses ('x').

In the next graph, we use two times fewer x values by changing the code line:

x = np.linspace(0, 4 * np.pi, 40)  # data: x values

to:

x = np.linspace(0, 4 * np.pi, 20)  # data: x values

This makes it possible to see visually how interpolation approximates the original function.

3.2.2 Least-Squares Approximation of a Solution

Least-squares approximation of a solution in the one-dimensional case is the process of fitting a linear model to data. In this example, we approximate a dataset by a straight line in the x-y plane.

y = wx + b

Here y is the value that we are trying to predict, x is the independent variable, w is the gradient, or slope, of the line, and b is the bias, or the point where the line intersects the Y axis.

We determine the gradient w and the bias b using data pairs $[x_i, y_i]$ . In this example, the equation can be rewritten as:

y = X \cdot w

In this equation, X is a data matrix that includes the values [x, 1], with the 1 acting as a constant that accounts for the bias b. Now w is a parameter vector containing the gradient and the bias. We can now use np.linalg.lstsq to solve for w:

import numpy as np  # import the NumPy package and assign it the alias np
import matplotlib.pyplot as plt

w_true, b_true = 2.0, 0.0  # set the true values for the gradient and bias
x = np.array([0, 1, 2, 3, 4, 5])  # create an array of x values
X = np.vstack([x, np.ones(len(x))]).T  # form matrix X

# generate y values with random noise added
y = w_true * x + b_true + np.random.rand(len(x)) / 0.1

w, b = np.linalg.lstsq(X, y, rcond=None)[0]  # compute parameters w and b
plt.plot(x, y, 'o', label='Original data', markersize=10)
plt.plot(x, w * x + b, 'r', label='Approximated line')
plt.legend()
plt.show()

The resulting model is the simplest form of machine learning, known as linear regression. It shows how dependencies between variables can be approximated by a linear function. This method is used in machine learning to predict values or to understand relationships between variables.

3.2.3 One-Dimensional Interpolation Using interp1d

One-dimensional interpolation is a simple but powerful tool. Let us consider the interp1d interpolation method from the SciPy library.

The interp1d function makes it possible to perform interpolation for one-dimensional data. This can be linear interpolation, where new points are created by connecting neighboring data points with straight lines, or a more complex form, for example cubic interpolation, where cubic polynomials are used to smoothly connect the points.

Consider the following example:

import numpy as np  # import the NumPy package and assign it the alias np
# Import the Matplotlib library for data visualization
import matplotlib.pyplot as plt
# Import the interp1d function for interpolation
from scipy.interpolate import interp1d

x0, xL = 0, 20  # set the initial (x0) and final (xL) values of the x range

# Generate 11 uniformly distributed points between x0 and xL
x = np.linspace(x0, xL, num=11, endpoint=True)

# Compute y values as the cosine of x cubed, scaled by 1/8
y = np.cos(-x ** 3 / 8.0)

print('Data x:', x)
print('Data y:', y)
print('x.shape:', x.shape, 'y.shape:', y.shape)  # print sizes of x and y arrays

f = interp1d(x, y)  # create a linear interpolation function
f2 = interp1d(x, y, kind='cubic')  # create a cubic interpolation function

# Generate new points for prediction
xnew = np.linspace(x0, xL, num=41, endpoint=True)

# Visualize the original data and the interpolation results
plt.plot(x, y, 'o', xnew, f(xnew), '-', xnew, f2(xnew), '-')
plt.legend(['data', 'linear', 'cubic'], loc='best')  # add legend
plt.show()  # show graph

Data x: [ 0.  2.  4.  6.  8. 10. 12. 14. 16. 18. 20.]
Data y: [ 1.          0.54030231 -0.14550003 -0.29213881  0.39185723  0.78771451
 -0.71798508 -0.84383778 -0.99683339  0.98869558  0.56237908]
x.shape: (11,) y.shape: (11,)

3.2.4 Representation of a Two-Dimensional Spline Using bisplrep

A two-dimensional spline is a means of creating smooth surfaces from a dataset in two dimensions, usually denoted as x and y. Each point in this two-dimensional space has a corresponding value, often denoted as z. Splines are used to create a continuous and smooth function that passes through these points or approximates their arrangement.

The bisplrep function, short for B-spline representation, is used to construct a B-spline, which is a mathematical model for describing smooth curves and surfaces. B-splines provide flexibility and accuracy when working with data, making it possible to control the shape of the approximation through a number of parameters.

Consider the following example:

# Import the required libraries
import numpy as np  # for working with arrays and mathematical functions
import matplotlib.pyplot as plt  # for data visualization
from scipy import interpolate  # for interpolation and working with splines

# Create a two-dimensional coordinate grid
# np.mgrid creates a two-dimensional grid in the range from -2 to 2,
# with 20 points along each axis
grid_x, grid_y = np.mgrid[-2:2:20j, -2:2:20j]

# Compute values according to the Gaussian function
# Apply the Gaussian function to each grid point
gaussian_z = np.exp(-0.5 * (grid_x ** 2 + grid_y ** 2))

# Visualize the original data
plt.figure()  # create a new figure
# draw a color map of the original Gaussian-function data
plt.pcolor(grid_x, grid_y, gaussian_z, shading='auto')
plt.colorbar()  # add a color bar to understand the values
plt.title("Gaussian function at discrete points")  # set graph title
plt.show()  # show graph

# Create a new, denser grid for interpolation
# Create a denser grid for detailed interpolation
fine_grid_x, fine_grid_y = np.mgrid[-2:2:100j, -2:2:100j]

# Create a B-spline using the bisplrep function
# bisplrep creates B-spline parameters from the original data
spline_tck = interpolate.bisplrep(grid_x, grid_y, gaussian_z, s=0)

# Compute spline values on the new grid
# bisplev computes spline values on the new grid
interpolated_z = interpolate.bisplev(fine_grid_x[:, 0], fine_grid_y[0, :], spline_tck)

# Visualize the interpolated data
plt.figure()  # create another figure
# draw a color map of the interpolated data
plt.pcolor(fine_grid_x, fine_grid_y, interpolated_z, shading='auto')
plt.colorbar()  # add a color bar
plt.title("Interpolated Gaussian function")  # set graph title
plt.show()  # show graph

The Gaussian function used in the code above is defined as follows:

def gaussian_function(x, y):
    return np.exp(-0.5 * (x ** 2 + y ** 2))

In this function, x and y are coordinates for which the value of the Gaussian function is computed. As a rule, the Gaussian function has a bell-shaped form centered at the point (0, 0) and decreases as the distance from the center increases.

In this specific case, the formula used represents a two-dimensional normal distribution without displacement and with identical standard deviations along both axes, centered at the point (0, 0). The value 0.5 in the exponent determines the “width” of the function’s bell; larger values lead to a “sharper” and higher bell, while smaller values lead to a “flatter” and wider one.

Experiment with the code; this will help you understand the material better and feel more confident when writing code.

3.2.5 Radial Basis Functions for Smoothing and Interpolation

To understand how radial basis functions (RBFs) work, let us examine their main aspects:

Basis functions for approximation: RBFs are used to approximate functions, that is, to find functions that approximately correspond to a dataset. This is especially useful when we have a set of data points and want to find a smooth function that passes through these points or close to them.

Distance functions: RBFs are based on the concept of distance. These functions depend on the distance from a central point, which makes them radially symmetric. This property allows RBFs to be very flexible when working with data distributed nonuniformly or in complex multidimensional spaces.

Resistance to overfitting: RBFs are usually less prone to overfitting than some other approximation methods. This is related to their smoothness and their ability to generalize.

The SciPy library provides different types of RBFs, each with its own features:

"multiquadric": sqrt((r/self.epsilon)**2 + 1)

This is a universal function that provides a smooth approximation. It can handle different types of data and is often used by default.

"inverse": 1.0/sqrt((r/self.epsilon)**2 + 1)

The inverse multiquadric function creates a smooth surface that is well suited for approximating smooth functions.

"gaussian": exp(-(r/self.epsilon)**2)

The Gaussian function provides a very smooth approximation and is often used in statistical analysis and machine learning.

"linear": r

The linear function is simple and may be useful for some basic approximation tasks.

"cubic": r**3

The cubic function provides a sharper approximation compared with the linear or Gaussian function.

"quintic": r**5

The quintic function provides an even higher degree of smoothness.

"thin plate": r*2 log(r)

The “thin plate” function is especially useful for interpolation of spatial data and is often used in geostatistics and computer graphics.

Many RBFs, for example the multiquadric or Gaussian functions, use an epsilon parameter that controls the “width” of the function. By changing its value, one can control the degree of smoothness or localization of the approximation. By selecting a suitable RBF type and tuning the parameters, one can achieve an optimal balance between approximation accuracy and resistance to overfitting.

Let us consider one-dimensional examples.

import numpy as np
from scipy.interpolate import Rbf, InterpolatedUnivariateSpline
import matplotlib.pyplot as plt

# Data generation
np.set_printoptions(formatter={'float': '{: 0.3f}'.format})

# create an array of points from -10 to 10, 11 points in total
sample_points = np.linspace(-10, 10, 11)

# compute the values as cosine of sample_points
sample_values = np.cos(sample_points)

# Create a finer grid for displaying interpolated data
# create a denser grid of points from -10 to 10
interpolation_points = np.linspace(-10, 10, 100)

# Use splines for interpolation
# create a spline interpolator
spline_interpolator = InterpolatedUnivariateSpline(sample_points, sample_values)

# compute interpolated values on the dense grid
spline_values = spline_interpolator(interpolation_points)

plt.subplot(2, 1, 1)  # prepare to display two plots together
# display original data as blue points
plt.plot(sample_points, sample_values, 'bo')
# display the original cosine function as a red line
plt.plot(interpolation_points, np.cos(interpolation_points), 'r')
# display spline-interpolated values as a green line
plt.plot(interpolation_points, spline_values, 'g')
plt.title('Interpolation using a one-dimensional spline')
# Add a legend to the graph
# plt.show()

# Use RBF for interpolation
# create an RBF interpolator
rbf_interpolator = Rbf(sample_points, sample_values)

# compute interpolated values using the RBF method
rbf_values = rbf_interpolator(interpolation_points)

plt.subplot(2, 1, 2)  # prepare the second plot
# display original data as blue points
plt.plot(sample_points, sample_values, 'bo')
# display the original cosine function as a red line
plt.plot(interpolation_points, np.cos(interpolation_points), 'r')
# display RBF-interpolated values as a green line
plt.plot(interpolation_points, rbf_values, 'g')
plt.title('Interpolation using RBF')
# Add a legend to the graph
plt.legend(['original data', 'original function', 'interpolated values'],
           loc='upper center', bbox_to_anchor=(0.5, -0.1),
           fancybox=True, shadow=True, ncol=4)
plt.show()

Now consider two-dimensional examples:

import numpy as np
from scipy.interpolate import Rbf
import matplotlib.pyplot as plt
from matplotlib import cm

# Generate 2D test data
num_points = 10  # number of points
x = np.random.rand(num_points) * 6 - 3  # random x in the range [-3, 3]
y = np.random.rand(num_points) * 6 - 3  # random y in the range [-3, 3]
z = np.sin(np.sqrt(x ** 2 + y ** 2))  # z as sine of distance from the origin

# Create a grid for interpolation
grid_size = 100
xi = np.linspace(-3, 3, grid_size)
yi = np.linspace(-3, 3, grid_size)
XI, YI = np.meshgrid(xi, yi)  # create grid

# Use RBF for interpolation
rbf_interpolator = Rbf(x, y, z, function='multiquadric', epsilon=0.5)  # create RBF interpolator
ZI = rbf_interpolator(XI, YI)  # interpolation

# Visualize the results
plt.figure(figsize=(8, 6))
plt.pcolor(XI, YI, ZI, cmap=cm.viridis, shading='auto')  # display interpolated data
plt.scatter(x, y, 30, z, cmap=cm.viridis, edgecolor='k')  # display original data points
plt.title('2D RBF interpolation - multiquadric')
plt.xlim(-3, 3)
plt.ylim(-3, 3)
plt.colorbar()
plt.show()

import numpy as np
from scipy.interpolate import Rbf
import matplotlib.pyplot as plt

# Generate random 3D data
x, y, z = np.random.rand(3, 30)  # 30 random points in each dimension in the range [0, 1]
d = np.sin(x) + np.cos(y) + np.tan(z)  # compute d values from a function of x, y, and z

# Create an RBF interpolator
rbf_interpolator = Rbf(x, y, z, d, function='multiquadric')

# Create a grid for interpolation
xi, yi, zi = np.meshgrid(np.linspace(0, 1, 10),
                         np.linspace(0, 1, 10),
                         np.linspace(0, 1, 10))

# Interpolation
di = rbf_interpolator(xi, yi, zi)

# Visualize original and interpolated data
fig = plt.figure(figsize=(12, 6))

# Visualize original data
ax1 = fig.add_subplot(121, projection='3d')
ax1.scatter(x, y, z, c=d, cmap='viridis')
ax1.set_title("Original data")
ax1.set_xlabel("X")
ax1.set_ylabel("Y")
ax1.set_zlabel("Z")

# Visualize interpolated data
ax2 = fig.add_subplot(122, projection='3d')
ax2.scatter(xi.flatten(), yi.flatten(), zi.flatten(), c=di.flatten(), cmap='viridis')
ax2.set_title("Interpolated data")
ax2.set_xlabel("X")
ax2.set_ylabel("Y")
ax2.set_zlabel("Z")

plt.show()

3.3 Singular Value Decomposition (SVD)

Singular value decomposition (SVD) is a matrix factorization method used to represent a matrix through orthogonal or unitary factors and singular values. In this section, SVD is considered as a tool for data analysis, dimensionality reduction, and compression in machine-learning tasks.

3.3.1 SVD Formulation

SVD is a matrix factorization method used in many fields, such as signal processing, statistics, and machine learning. The essence of SVD is to decompose a matrix into three other matrices; it is applicable to both square and rectangular matrices.

The SVD decomposition of a matrix A of size m × p is represented as:

A = U\Sigma V^{*}

where:

U is a unitary matrix of size m × m.

Σ (sigma) is a rectangular diagonal matrix of size m × p with nonnegative numbers on the diagonal. These numbers are called the singular values of matrix A.

V* is the conjugate-transpose unitary matrix of size p × p.

The unitary matrices U and V have the property of preserving length and angle, which makes them ideal for describing rotations and reflections in space. In the context of SVD, these matrices provide bases for the input and output spaces of matrix A.

The singular values Σ indicate the directions in which matrix A “stretches” space most strongly. For example, if one singular value is much larger than the others, this indicates that most of the “energy” or “information” of the matrix is concentrated along the corresponding singular vector.

The left singular vectors, the columns of U, represent directions in the original space that, after transformation by matrix A, turn into the corresponding right singular vectors, the columns of V. These vectors form orthonormal bases in their respective spaces.

SVD is a powerful tool in numerical analysis that allows us to “disassemble” a matrix into its component parts and understand how different parts of this matrix interact with one another.

3.3.2 Algorithms for SVD

One frequently used method for computing SVD is based on QR decomposition. QR decomposition is a process in which a matrix is decomposed into an orthogonal matrix Q and an upper triangular matrix R. This approach is well suited for working with square and non-square matrices.

Another approach to computing SVD is based on eigenvalue decomposition. This method is theoretically simple, but in practice it is rarely used for large matrices because of numerical-stability problems. For our theoretical analysis and formula derivation, however, this is not a problem, and therefore we consider this approach. To make this simple approach work, we need to impose several conditions on matrix A.

Suppose we have a matrix A of size m × p, consisting of real numbers, with m > p and rank equal to p.

Formation of the normal matrix B: this matrix will be symmetric and square, of size p × p.

B = A^T A

Then matrix B is subjected to eigenvalue decomposition, leading to the formula:

B = V_e \Lambda V_e^T

where:

V_e is an orthonormal matrix of eigenvectors,

Λ is a diagonal matrix of eigenvalues.

It is known that A has an SVD decomposition, that is:

A = U\Sigma V^T

Because matrix A has rank p, all singular values in Σ are positive real numbers. Using the formula $A = U\Sigma V^T$ , we have:

A^T A = (U\Sigma V^T)^T(U\Sigma V^T) = V\Sigma U^T U\Sigma V^T = V\Sigma^2 V^T = B

Since $U^T U = I$ , and since Σ is diagonal and does not change under transposition, comparing the previous formula with $B = V_e \Lambda V_e^T$ gives:

V = V_e

\Sigma = \sqrt{\Lambda}

Now using $A = U\Sigma V^T$ and the orthonormal property of V, namely $V^T V = I$ , we have:

AV = U\Sigma

Since Σ has all positive eigenvalues, we finally obtain:

U = AV\Sigma^{-1}

The statement that U is unitary follows from the fact that A is represented by a singular decomposition, where U and V are unitary matrices. Finally, if A is underdetermined, that is, if rank(A) < p, then matrix B will have zero eigenvalues. In such cases, we simply discard all zero eigenvalues and their corresponding eigenvectors. This still gives us SVD in reduced form, and the entire process described above remains valid.

The main problem with this method is numerical stability, especially for large matrices. The problem arises because condition numbers are squared during the formation of matrix $B = A^T A$ .

Although this method is useful for theoretical understanding of SVD and for proving that every matrix has an SVD, it is rarely used in practice for computations. In practice, other methods, such as QR decomposition, are more often used to compute SVD; these avoid the problems associated with forming matrix B.

3.3.3 Example

Let us consider an example of using SVD, or singular value decomposition, with Python and the NumPy library.

import numpy as np  # import the NumPy library

# Create matrix 'my_matrix' with random integer values of size 4x5
my_matrix = np.random.randint(0, 100, size=(4, 5))
print('-' * 20)
print("Original matrix:\n", my_matrix)
print('-' * 20)

# Perform singular value decomposition (SVD) of matrix 'my_matrix'
# U - left singular vectors,
# sigma - singular values,
# Vt - right singular vectors (transposed)
U, sigma, Vt = np.linalg.svd(my_matrix, full_matrices=False)

# Print the dimensions of U, sigma, and Vt for checking
print("\nShapes: U, sigma, Vt:", U.shape, sigma.shape, Vt.shape)
print('-' * 20)

# Print matrices U, sigma, and Vt
print("\nleft singular vectors U:\n", U)
print("\nsingular values Sigma:", sigma)
print("\nright singular vectors Vt:\n", Vt)
print('-' * 20)

# Convert the singular-value vector 'sigma' into a diagonal matrix 'sigma_mat'
sigma_mat = np.diag(sigma)
print("\ndiagonal matrix sigma:\n", sigma_mat)
print('-' * 20)

# Check whether the original matrix can be reconstructed
reconstructed_matrix = np.dot(U, np.dot(sigma_mat, Vt))
print("\nHas the original matrix been reconstructed? ",
      np.allclose(my_matrix, reconstructed_matrix))
print('-' * 20)

# Print the original and reconstructed matrices for comparison
print("\noriginal matrix:\n", my_matrix)
print("\nreconstructed matrix:\n", reconstructed_matrix)

--------------------
Original matrix:
 [[58 83 22 94 36]
 [22 54 45 82 91]
 [10 34 62 75 24]
 [55 26 46 52 68]]
--------------------

Shapes: U, sigma, Vt: (4, 4) (4,) (4, 5)
--------------------

left singular vectors U:
[[ 0.56147956  0.79253797  0.237659   -0.01193577]
 [ 0.57205124 -0.43106941  0.05107003 -0.69593706]
 [ 0.40034304 -0.01747367 -0.87370638  0.27578483]
 [ 0.44409856 -0.43099582  0.42136351  0.66292675]]

singular values Sigma: [242.50494188  59.13289951  48.44113155  33.06221611]

right singular vectors Vt:
[[ 0.30341564  0.42329775  0.34368231  0.63011555  0.46216321]
 [ 0.21314985  0.51921871 -0.38678068  0.26091348 -0.68359633]
 [ 0.60580112  0.07706081 -0.56275368 -0.35278672  0.43118031]
 [ 0.70218915 -0.36169658  0.4843455  -0.09172851 -0.36483051]]
--------------------

diagonal matrix sigma:
[[242.50494188   0.           0.           0.        ]
 [  0.          59.13289951   0.           0.        ]
 [  0.           0.          48.44113155   0.        ]
 [  0.           0.           0.          33.06221611]]
--------------------

Has the original matrix been reconstructed? True
--------------------

original matrix:
[[58 83 22 94 36]
 [22 54 45 82 91]
 [10 34 62 75 24]
 [55 26 46 52 68]]

reconstructed matrix:
[[58. 83. 22. 94. 36.]
 [22. 54. 45. 82. 91.]
 [10. 34. 62. 75. 24.]
 [55. 26. 46. 52. 68.]]

This code example and the comments should help clarify how SVD works and how it can be applied in real tasks.

3.3.4 SVD for Data Compression

Image compression using SVD: for this example we use an image loaded from the PIL, or Python Imaging Library, and perform SVD on it for data compression.

import numpy as np
from PIL import Image
import matplotlib.pyplot as plt

# Load image
# Replace 'path_to_image.jpg' with the path to your image
original_image = Image.open('path_to_image.jpg')

# Convert the image to black-and-white format to simplify processing
gray_image = original_image.convert('L')

# Convert the PIL image to a NumPy array for further processing
gray_array = np.array(gray_image)

print('resulting data array from the image\n', gray_array)
print('-' * 20)
print("Shape of the original image array:", gray_array.shape)

# Create a window to display images
plt.figure(figsize=(12, 6))

# Display the original image
plt.subplot(1, 2, 1)  # 1 row, 2 columns, position 1
plt.imshow(gray_array, cmap='gray')
plt.title("Original image")

# Perform SVD on the image matrix
U, S, Vt = np.linalg.svd(gray_array, full_matrices=False)
print('-' * 20)

# Print the number of singular values in the original image
print("Number of singular values in the original image:", len(S))

# Determine the number of singular values to use
k = 112  # experiment with this value for different compression levels
print('-' * 20)
print("Number of singular values used:", k)

# Build an approximate image
# Matrix multiplication to reconstruct the image using k singular values
approximated_image = np.dot(U[:, :k], np.dot(np.diag(S[:k]), Vt[:k, :]))

# Display the reconstructed image
plt.subplot(1, 2, 2)  # 1 row, 2 columns, position 2
plt.imshow(approximated_image, cmap='gray')
plt.title(f"Compressed image with {k} singular values")
plt.show()

resulting data array from the image
[[40 40 40 ... 74 74 73]
 [40 40 40 ... 74 74 74]
 [41 41 41 ... 74 74 75]
 ...
 [53 53 53 ... 61 61 61]
 [53 53 53 ... 61 61 61]
 [53 53 53 ... 61 61 61]]
--------------------
Shape of the original image array: (1125, 2000)
--------------------
Number of singular values in the original image: 1125
--------------------
Number of singular values used: 112

The variable k sets the number of singular values that will be used for approximate reconstruction of the image. This is a key parameter that affects image quality and the degree of compression. The expression np.dot(U[:, :k], np.dot(np.diag(S[:k]), Vt[:k, :])) multiplies the matrices U, S, and $V^T$ to reconstruct the image using only k singular values. Thus, we obtain a compressed version of the original image.

This example demonstrates how SVD can be used for image compression, reducing the amount of data needed for storage and transmission. It is important to note that the value of k, the number of singular values used, is a compromise between image quality and compression ratio. Too small a value of k can lead to significant loss of detail, while too large a value results in insignificant compression. By experimenting with different values of k, one can find the optimal balance for a particular task.

3.4 Principal Component Analysis

Principal Component Analysis (PCA) is a statistical method used in machine learning to reduce data dimensionality while preserving as much information as possible. First proposed by Karl Pearson in 1901, PCA makes it possible to transform a large set of variables that are often correlated with one another into a smaller set of uncorrelated variables called principal components.

PCA starts with an original dataset consisting of p variables, or features. The goal is to transform these data so that the new variables, the principal components, are linearly uncorrelated. These components are ordered so that the first principal component explains the largest share of the total variability in the data, the second explains the next most important part, and so on.

The transformation used in PCA is orthogonal. This means that the principal components are orthogonal, that is, uncorrelated, to one another in the multidimensional feature space. Orthogonal transformation ensures that each principal component represents a unique aspect of the data without information redundancy.

Principal components are ranked by the degree of explained variability in the data. The first principal component has the highest variability, the second has less, and so on. This ranking makes it possible to determine how many components should be retained to represent the data adequately while reducing the dimensionality of the original feature space.

One of the main applications of PCA in machine learning is data-dimensionality reduction. This is especially useful when working with datasets that have a large number of features, because dimensionality reduction can help speed up model training and reduce the risk of overfitting.

When using PCA, it is important to remember that it is sensitive to variable scaling. Different feature scales can significantly affect analysis results, so it is often recommended to standardize the data before applying PCA. In addition, PCA relies on linear relationships between variables and may be ineffective for identifying nonlinear structures in data.

3.4.1 PCA Formulation

Consider a general matrix A of size m × p, consisting of real numbers, where m > p. First form matrix B:

B = A^T A

This leads to the creation of a symmetric square matrix of size p × p. This matrix is at least positive semidefinite and often positive definite (SPD). Next, eigenvalue decomposition of matrix B is performed:

B = V\Sigma V^T

where:

V is an orthonormal matrix of size p × p, consisting of p eigenvectors of matrix B.

Σ is a diagonal matrix of size p × p, containing eigenvalues that are nonnegative real numbers.

Then PCA is defined as:

A_{PCA} = AV

Thus, we project the original data A onto orthonormal eigenvectors. The resulting matrix $A_{PCA}$ preserves the shape of the original matrix A, namely m × p.

Using all eigenvectors, the original matrix A can be reconstructed:

A_r = A_{PCA}V^T = AVV^T = A

This is possible because the eigenvectors are orthonormal. However, often only the first few eigenvectors, ranked by eigenvalues, are used, because they contain most of the information of the original matrix A. An example of partial reconstruction using k eigenvectors is:

A_r = A_{PCA}[0:m,0:k]V^T[0:k,0:p] = A[0:m,0:p]V[0:p,0:k]V^T[0:k,0:p] \neq A

In general, this is not equal to the original A, but it can often be very close to it. In this case, the storage size is m × k + k × p, which may be much smaller than the original size m × p.

If the dimensions of matrix A are such that m < p, its transpose can be used and an analogous process applied. In this case, matrix B is formed as:

B = A A^T

This leads to the creation of a symmetric square matrix of size m × m. After eigenvalue decomposition of this matrix is performed, we obtain:

B = V\Sigma V^T

where V and Σ are defined analogously to the previous case, but with dimensions m × m.

Then PCA is defined as:

A_{PCA} = V^T A

The result preserves the shape of the original matrix A, namely m × p. Full reconstruction of the original data is also possible using all eigenvectors:

A_r = VA_{PCA} = VV^T A = A

We can use only a small number of eigenvectors to reconstruct matrix A. For example, if we use k ≤ m eigenvectors, we obtain:

A_r = V[0:k,0:m]V^T[0:m,0:k]A_{PCA}[0:k,0:p] \neq A

For large systems, instead of direct formation of matrix B and eigenvalue decomposition, more stable algorithms are used, such as QR transformation, which was discussed above.

3.4.2 Examples

Example 1:

Let us consider an example of PCA used to analyze a dataset of height and weight measurements. This example helps illustrate how PCA can be used to identify the main directions of variation in data.

Imagine that we have a dataset containing height measurements in centimeters and weight measurements in kilograms for a group of people. Our goal is to use PCA to identify the main directions of variation in these two variables.

import numpy as np
import matplotlib.pyplot as plt

# Original data
height = np.array([170, 168, 177, 181, 172, 171, 169, 175, 174, 178])
weight = np.array([70, 66, 78, 85, 72, 68, 65, 74, 71, 80])

# Create the data matrix
A = np.vstack((height, weight)).T

PCA begins with centering the data, that is, subtracting the mean value of each feature from the corresponding values.

# Center the data
mean_height = np.mean(height)
mean_weight = np.mean(weight)
A_centered = A - np.array([mean_height, mean_weight])

The next step is to compute the covariance matrix, which shows how the variables change together.

# Compute the covariance matrix
cov_matrix = np.cov(A_centered.T)

Now we can perform PCA by finding the eigenvalues and eigenvectors of the covariance matrix.

# Perform PCA
eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)

# Sort eigenvalues and eigenvectors
index = eigenvalues.argsort()[::-1]
eigenvalues = eigenvalues[index]

Eigenvalues show how much variation is explained by each principal component. Eigenvectors show the directions of the principal components in the original feature space.

# First principal component
first_component = eigenvectors[:, 0]

# Visualization
plt.scatter(A_centered[:, 0], A_centered[:, 1])
plt.quiver(mean_height, mean_weight, first_component[0], first_component[1], scale=5, color='r')
plt.xlabel('Height')
plt.ylabel('Weight')
plt.title('PCA of height and weight')
plt.axis('equal')
plt.show()

In this example, the first principal component, shown by the red arrow, indicates the direction of greatest variation in the data. If height and weight are strongly correlated, this component will point in the direction in which both measurements increase together. This may represent, for example, the general tendency for weight to increase with height.

Now that we have determined the principal components, we can project the original data onto these components for dimensionality reduction or further analysis.

# Project data onto the first principal component
projected_data = np.dot(A_centered, eigenvectors)

# Visualize the projection
plt.scatter(projected_data[:, 0], projected_data[:, 1])
plt.xlabel('First principal component')
plt.ylabel('Second principal component')
plt.title('Projection of height and weight data onto principal components')
plt.axis('equal')
plt.show()

In this graph, the data points are now represented in a coordinate system defined by the principal components. The first axis, the x axis, corresponds to the first principal component and captures the greatest variation in the data. The second axis, the y axis, is the second principal component, which captures the remaining variation.

This example illustrates how PCA can be used to identify the main directions of variation in data. PCA is especially useful in situations where there are many variables and it is necessary to find the most significant directions of variation to simplify analysis or prepare data for machine learning.

Example 2

Let us write code for truncated PCA, or Principal Component Analysis, which can be used for image compression. In this example, we load an image, convert it to grayscale, apply PCA, and then reconstruct the image using a limited number of principal components.

import numpy as np
import matplotlib.pyplot as plt
from PIL import Image

# Define a function to perform PCA
def princomp(A, numpc=0):
    # Subtract the mean value for each feature to center the data
    A = A - np.mean(A, axis=0)

    # Compute the covariance matrix of the data
    cov_matrix = np.dot(A.T, A) / (A.shape[0] - 1)

    # Find eigenvalues and eigenvectors of the covariance matrix
    eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)

    # Convert eigenvalues and eigenvectors to real numbers
    eigenvalues = np.real(eigenvalues)
    eigenvectors = np.real(eigenvectors)

    # Sort eigenvalues and corresponding eigenvectors in descending order
    idx = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]

    # Truncate the number of principal components to the specified number numpc
    if numpc < A.shape[1] and numpc >= 0:
        eigenvectors = eigenvectors[:, :numpc]

    # Project the data onto principal components to obtain new features
    score = np.dot(A, eigenvectors)
    return eigenvectors, score, eigenvalues

# Load image from the specified path
image = Image.open('path_to_image.jpg')  # replace with the path to the image

# Convert the image to grayscale
image = image.convert('L')

# Convert the image to a two-dimensional NumPy array
image_array = np.array(image)

# Apply PCA to the image, limiting the number of principal components to 50
eigenvectors, score, eigenvalues = princomp(image_array, numpc=50)

# Reconstruct the image from the compressed representation
reconstructed_image = np.dot(score, eigenvectors.T) + np.mean(image_array, axis=0)

# Use only the real part of the data for the reconstructed image
reconstructed_image = np.real(reconstructed_image)

# Create a new figure to display the images
plt.figure(figsize=(8, 4))

# Display the original image
plt.subplot(1, 2, 1)
plt.imshow(image_array, cmap='gray')  # show the original image
plt.title('Original image')  # title for the original image
plt.axis('off')  # turn off axes for the original image

plt.subplot(1, 2, 2)
# Show the reconstructed image in grayscale
plt.imshow(reconstructed_image, cmap='gray')
plt.title('Reconstructed image')  # title for the reconstructed image
plt.axis('off')  # turn off axes for the reconstructed image
plt.show()

In this code, we first define the princomp function, which performs PCA on the input matrix A. Then we load an image, convert it to grayscale, and apply PCA to it. As a result of PCA, we obtain eigenvectors, data projections onto the principal components, called score, and eigenvalues. At the end, we reconstruct the image from its compressed representation and visualize the original and reconstructed images.

We see that using 50 principal components provides fairly good image quality compared with the original. But if we use only 5 principal components, we obtain a readable image with lower quality.

As an important stage in learning Principal Component Analysis (PCA), it is strongly recommended to carry out independent experiments with the code. This not only deepens theoretical understanding of the method, but also develops practical skills in its application. By testing different code variants and adapting it to various tasks, you will gain valuable experience that contributes to a deeper and more intuitive understanding and helps you navigate real data-analysis scenarios. Experimenting with code also makes it possible to visualize and interpret the results.

3.5 Finding Numerical Roots

In machine learning, especially in optimization algorithms, it is often necessary to find roots of nonlinear equations. The scipy.optimize module in Python provides the fsolve() function, intended for finding roots of a set of nonlinear equations defined as f(x) = 0. This function is especially useful when the expected locations of the roots are known.

The fsolve() function is a wrapper for algorithms in MINPACK, which use Newton’s iteration method. This method consists of choosing an initial approximation to the root and then refining it by local linearization of the function.

The application of Newton’s method can be demonstrated visually by an example. Consider a function of one variable:

from scipy.optimize import fsolve

def stress_strain_curve(strain):
    # Example function
    return strain ** 3 - 0.2 * strain ** 2 + 0.1 * strain - 0.05

# Find the root of the function
strain_root = fsolve(stress_strain_curve, 0.1)
print('Root of the equation:', strain_root)

Root of the equation: [0.33940705]

Similarly, fsolve() can be applied to a system of equations with several variables, which is often encountered in multidimensional problems:

from scipy.optimize import fsolve

def market_equilibrium(x):
    # Two functions representing market demand and supply
    demand = 10 - x[0] + 0.5 * x[1]
    supply = x[0] - 5 - 0.25 * x[1]
    return [demand, supply]

# Initial guess for price and quantity
initial_guess = [2, 2]
price_quantity = fsolve(market_equilibrium, initial_guess)
print('Market equilibrium:', price_quantity)

Market equilibrium: [-4.03896783e-27 -2.00000000e+01]

It is important to emphasize that in machine learning, when working with polynomials or algebraic equations, complex roots may appear even when the coefficients of the equations are real numbers. This is because complex space is geometrically closed, unlike real space. Thus, a polynomial of degree n must have n roots, which may lie either in real space or in complex space.

3.6 Numerical Integration in the Context of Machine Learning

Numerical integration considers the task of estimating a definite integral of a function. The main idea is to approximate the integral by a sum of function values at certain points multiplied by small intervals.

The trapezoidal method: this is one of the basic methods used in numerical integration. The method consists of dividing the area under the graph of a function into small trapezoids and calculating their total area. This is a simple and effective method for functions that do not change their behavior too quickly.

Simpson’s method: this method is based on approximating a function by parabolas on small intervals and computing the area under these parabolas. Simpson’s method provides higher accuracy than the trapezoidal method, especially for functions that change in a more complex way.

The NumPy library provides convenient tools for implementing these methods. For example, the numpy.trapz() function can be used to apply the trapezoidal method.

Example of use:

import numpy as np

# Function that must be integrated
def function(x):
    return np.sin(x)

# Generate points and function values
x = np.linspace(0, np.pi, 100)
y = function(x)

# Apply the trapezoidal method
integral = np.trapz(y, x)
print("Integral value:", integral)

Integral value: 1.9998321638939927

3.6.1 Trapezoidal Rule

Let us consider the trapezoidal rule in more detail. It is a method of numerical integration used for approximate computation of the definite integral of a function. The formula for the trapezoidal rule is:

\int_a^b f(x)\,dx \approx \frac{1}{2}\sum_{k=1}^{n_s}(x_k-x_{k-1})(f(x_k)+f(x_{k-1}))

To demonstrate the trapezoidal rule in practice, we define a simple polynomial function and sample its values in the finite range [a, b] at $n_s$ uniformly distributed points.

Below is Python code demonstrating the use of the trapezoidal rule:

# Import the required libraries
import numpy as np  # import the NumPy library for numerical computations
import matplotlib.pyplot as plt  # import the Matplotlib library
from scipy.integrate import quad  # import the quad function from SciPy

# Configure output format
# Configure the NumPy output format for floating-point numbers
np.set_printoptions(formatter={'float': '{: 0.3f}'.format})

# Define a polynomial function
def f(x):  # define function f(x)
    return 5 * x ** 3 - 6 * x ** 2 - 7 * x + 6

# Sample function values
# Set the starting point (a), ending point (b), and number of points (n)
a, b, n = -1., 2, 400

# Create an array of n uniformly distributed points between a and b
x = np.linspace(a, b, n)
y = f(x)  # compute f(x) at each point of array x

# Function integration
ns = 6  # set the number of points for the trapezoidal method
# Create an array of ns uniformly distributed points for integration
xint = np.linspace(a, b, ns)
yint = f(xint)  # compute f(x) at each point of array xint

# Visualize the results
plt.plot(x, y, lw=2)  # plot function f(x)
# Fill the area under the graph between xint and yint
plt.fill_between(xint, 0, yint, facecolor='gray', alpha=0.4)
# Add text with the integral to the plot
plt.text((a + b) / 2, 12, r"$\int_a^b f(x)dx$",
         horizontalalignment='center', fontsize=15)
plt.show()  # display graph

# Compute the integral
# Compute the exact integral value and the error
integral, error = quad(f, a, b)
# Compute the approximate integral value by the trapezoidal method
integral_trapezoid = sum((xint[1:] - xint[:-1]) * (yint[1:] + yint[:-1])) / 2

# Print results
# Print the exact integral value and its error
print("Exact result:", integral, "+/-", error)
# Print the approximate integral value obtained by the trapezoidal method
# and the number of points used
print("Trapezoidal approximation result with", len(xint), "points:", integral_trapezoid)

Exact result: 8.25 +/- 1.2419279516758014e-13
Trapezoidal approximation result with 6 points: 8.520000000000001

This example shows that the results obtained with the trapezoidal rule are approximate. The computation error may be significant, especially with a small number of sample points. The accuracy of trapezoidal-rule results directly depends on the number of intervals used, $n_s$ . The larger the number of intervals, the more accurate the approximation; however, increasing the number of intervals also increases computational complexity.

The trapezoidal rule is a simple and effective method of numerical integration. This method makes it possible to obtain approximate integral values, which is often sufficient for practical purposes in machine learning. However, it is always important to remember the potential approximation error and to seek a balance between computational accuracy and computational complexity.

3.6.2 Gaussian Quadrature

Gaussian quadrature, also known as Gaussian integration, is a highly efficient technique for numerical integration. This method is based on selecting and summing function values at specially chosen points, called Gaussian points, using the corresponding weights. The approach provides high computational accuracy, especially when working with polynomial integrands, because Gaussian points correspond to the roots of Legendre polynomials on the interval [−1, 1].

Gaussian integration is used when computing integrals in methods based on Bayesian statistics or when solving optimization problems.

Consider an example in which we use Gaussian quadrature to compute the integral of a function. For this example, we take a function that may represent, for example, the probability of classifier error.

import numpy as np  # import the NumPy library
from scipy.integrate import quad  # import the quad function from SciPy

def error_function(x):
    """Representation of a classifier error function"""
    # Returns e raised to the power (-x^2), which is an exponential function
    return np.exp(-x ** 2)

# Set the integration interval
# The values a and b define the lower and upper limits of integration
a, b = 0, 1

# Integrate the function over the specified interval
# The quad function is used for numerical integration
integral, error = quad(error_function, a, b)
# error_function from a to b.
# Returns the integral value and the estimated error

# Print the integral value
print(f"Integral of the classifier error function on interval [{a}, {b}]:", integral)
# Print the estimate of the integration error
print("Estimated computational error:", error)

Integral of the classifier error function on interval [0, 1]: 0.7468241328124271
Estimated computational error: 8.291413475940725e-15

For general complex integrand functions, Gaussian quadrature may not give an exact solution. Its accuracy, however, will still be much better than that of the trapezoidal rule or the rectangle rule, which we have not discussed but which is very similar to the trapezoidal rule. In other words, to obtain solutions of similar accuracy, Gaussian integration uses fewer sampling points.

This approach to integration can be used to estimate the error probability in classification algorithms or to integrate functions that arise in Bayesian methods. The efficiency of Gaussian quadrature makes it possible to reduce computational costs, which is especially valuable when working with large amounts of data or complex models.

3.7 Data Preprocessing

Data preprocessing is a key stage in preparation for analysis and training machine-learning models. In this section, we consider the main methods and techniques applied to data before model training begins.

Data are often represented as a dataset X, where $X \in \mathbb{X}^{m \times p}$ . Here m denotes the number of data points, or observations, and p denotes the number of features or variables. In real problems, data values may vary greatly, which can lead to numerical-stability problems during model training.

To ensure stability and improve the performance of machine-learning models, data are usually subjected to preprocessing that includes normalization or scaling. There are two main methods of data normalization: min-max scaling and standard scaling.

3.7.1 Min-Max Scaling

Min-max scaling is a preprocessing stage in machine learning. This method standardizes the range of features, which is especially important for algorithms sensitive to feature scale, such as gradient descent.

The scaling formula is:

X_{scaled} = \frac{X - X.\min(axis=0)}{X.\max(axis=0) - X.\min(axis=0)}

where $X.\min$ and $X.\max$ are the minimum and maximum feature values. This method transforms all values of each feature into the range from 0 to 1. The method can be generalized to transform values into an arbitrary range [a, b] as follows:

X_{scaled} = a + \frac{X - X.\min(axis=0)}{X.\max(axis=0) - X.\min(axis=0)}

This makes it possible to adapt scaling to specific needs.

Once this scaling transformation has been performed on the training dataset, $X.\min$ and $X.\max$ can be used to perform exactly the same transformation on the test dataset to ensure consistency for correct predictions.

Below is simple code for performing min-max scaling.

import numpy as np

# Set output parameters for better readability
np.set_printoptions(precision=4)

# Modified dataset for the example
data_set = [
    [4, -1, 2],
    [9, 1, 5],
    [14, 2, 8],
    [19, 4, 11]
]

print('-' * 20)
print(f"Original training dataset:\n{data_set}")
print('-' * 20)

# Convert the list to a NumPy array for convenient computations
data_array = np.array(data_set)

# Perform min-max scaling
data_min = data_array.min(axis=0)
data_max = data_array.max(axis=0)
data_normalized = (data_array - data_min) / (data_max - data_min)

print(f"Scaled training dataset:\n{data_normalized}")
print(f"Maximum values for each feature:\n{data_max}")
print(f"Minimum values for each feature:\n{data_min}")
print('-' * 20)

# Modified test dataset
test_data = np.array([[0, 4, 6], [7, 5, 6]])
test_data_normalized = (test_data - data_min) / (data_max - data_min)
print("Scaled test dataset:\n", test_data_normalized)
print('-' * 20)

# Inverse transformation for the training dataset
data_original = data_normalized * (data_max - data_min) + data_min
print("Training dataset transformed back:\n", data_original)
print('-' * 20)

# Inverse transformation for the test dataset
test_data_original = test_data_normalized * (data_max - data_min) + data_min
print("Test dataset transformed back:\n", test_data_original)

--------------------
Original training dataset:
[[4, -1, 2], [9, 1, 5], [14, 2, 8], [19, 4, 11]]
--------------------
Scaled training dataset:
[[0.     0.     0.    ]
 [0.3333 0.4    0.3333]
 [0.6667 0.6    0.6667]
 [1.     1.     1.    ]]
Maximum values for each feature:
[19  4 11]
Minimum values for each feature:
[ 4 -1  2]
--------------------
Scaled test dataset:
[[-0.2667  1.      0.4444]
 [ 0.2     1.2     0.4444]]
--------------------
Training dataset transformed back:
[[ 4. -1.  2.]
 [ 9.  1.  5.]
 [14.  2.  8.]
 [19.  4. 11.]]
--------------------
Test dataset transformed back:
[[0. 4. 6.]
 [7. 5. 6.]]

It is clearly visible that min-max scaling does not damage the dataset. All data can be transformed back if necessary.

3.7.2 One-Hot Encoding

Many machine-learning datasets use categorical features. For example, the variable “color” may take the values “red,” “green,” and “blue.” To build a machine-learning model, it is necessary to transform these categorical values into numerical ones. Consider a one-dimensional feature vector initially given as [[green], [red], [0], [blue]]. One encoding method is to assign each color a unique number, for example [[1], [2], [0], [3]].

However, such direct encoding may implicitly introduce an order relationship between categories, which is not always desirable. For example, a model may interpret “blue” (3) as “greater” or “more important” than “green” (1), which may be incorrect.

To solve this problem, the one-hot encoding method is often used. In this approach, each category is transformed into a separate feature with a binary value. Scaling this type of data usually does not change it, because the values are already represented in the 0 and 1 format. The inverse transformation makes it possible to restore the original values.

One-hot encoding is an effective solution for working with categorical data. It avoids the mistaken interpretation of categorical variables as numerical variables with order relationships, which can lead to incorrect conclusions during data analysis. It is important to choose appropriate data-preprocessing methods depending on the nature of the data and the requirements of the task. One-hot encoding increases the dimensionality of the feature space, but it makes it possible to process data more accurately and correctly, ensuring that each category has equal influence on the model. This is especially important in algorithms that depend on distances between data points, such as k-nearest neighbors or clustering.

3.7.3 Standard Scaling

When the distribution of a dataset is close to a normal distribution, standard scaling can be used. The formula for standard scaling is as follows:

X_{scaled} = \frac{X - X.mean(axis=0)}{X.std(axis=0)}

Below is simple code for performing standard scaling.

import numpy as np

X = [[3, 5, 8],  # assumed training dataset
     [2, 6, 4.5],
     [3, 6, -6],
     [11, 5, 5]]

print('-' * 20)
print(f"Original training dataset X:\n{X}")
print('-' * 20)

X = np.array(X)
X_scaled = (X - X.mean(axis=0)) / X.std(axis=0)

print(f"Standard-scaled training dataset X:\n{X_scaled}")
print(f"Mean value for each feature:\n{X.mean(axis=0)}")
print(f"Standard deviation for each feature:\n{X.std(axis=0)}")
print('-' * 20)

Xtest = [[-15, 9, 2],  # assumed test dataset
         [4, 4, 4]]

Xt_scaled = (Xtest - X.mean(axis=0)) / X.std(axis=0)
print(f"Standard-scaled test dataset Xtest:\n{Xt_scaled}")
print('-' * 20)

X_back = X_scaled * X.std(axis=0) + X.mean(axis=0)
print(f"Training dataset transformed back:\n{X_back}")
print('-' * 20)

Xt_back = Xt_scaled * X.std(axis=0) + X.mean(axis=0)
print(f"Test dataset Xtest transformed back:\n{Xt_back}")

--------------------
Original training dataset X:
[[3, 5, 8], [2, 6, 4.5], [3, 6, -6], [11, 5, 5]]
--------------------
Standard-scaled training dataset X:
[[-0.48189987 -1.          0.96772426]
 [-0.75727123  1.          0.3068394 ]
 [-0.48189987  1.         -1.67581519]
 [ 1.72107098 -1.          0.40125152]]
Mean value for each feature:
[4.75  5.5   2.875]
Standard deviation for each feature:
[3.63145976 0.5        5.29593004]
--------------------
Standard-scaled test dataset Xtest:
[[-5.43858429  7.         -0.16522122]
 [-0.20652852 -3.          0.21242728]]
--------------------
Training dataset transformed back:
[[ 3.   5.   8. ]
 [ 2.   6.   4.5]
 [ 3.   6.  -6. ]
 [11.   5.   5. ]]
--------------------
Test dataset Xtest transformed back:
[[-15.   9.   2.]
 [  4.   4.   4.]]

It is important to take into account that scaling applied to features in the training dataset can also be applied to labels or target variables. This is especially relevant when the labels are numerical values rather than categories or a probability distribution. For example, if labels in a training dataset are numerical values such as weight or price, they can be scaled to improve model-training performance and stability.

However, it is important to remember that when testing a model or using it for prediction, labels should be returned to their original scale. This ensures the correct interpretation of results, because scaled labels may be unclear or misleading.

Check yourself

Which idea best describes the focus of "Basic Mathematical Computations"?

In machine learning, theoretical definitions are useful to check with numerical examples and visualizations.

Which actions help reinforce the chapter material?

Take quiz