Unit 5: Matrices - Subjective Questions

Random number generation is used in simulations, testing, games, and numerical experiments.

• Python's random module provides functions such as random(), randint(), uniform(), and choice().
• NumPy provides efficient functions through numpy.random, such as rand(), randn(), randint(), and uniform().
• A seed initializes the pseudorandom number generator. Using the same seed produces the same sequence, which makes experiments reproducible.

Example:

import random
import numpy as np

random.seed(10)
print(random.randint(1, 10))

rng = np.random.default_rng(10)
print(rng.integers(1, 11, size=5))

The generated values are called pseudorandom because they are produced by deterministic algorithms rather than by a truly random physical process.

NumPy is a Python library used for numerical computing. Its central data structure is the multidimensional homogeneous array called ndarray.

Important characteristics include:

• All elements generally have the same data type.
• Arrays can have one or more dimensions.
• Operations are implemented efficiently in compiled code.
• Mathematical operations can be applied to complete arrays without explicit Python loops.
• NumPy arrays use memory more efficiently than ordinary Python lists.

For example:

import numpy as np

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
print(a + b)  # [5 7 9]

In contrast, adding two Python lists usually concatenates them rather than performing element-wise addition. NumPy is therefore especially useful for matrices, scientific computing, statistics, and machine learning.

NumPy uses zero-based indexing, so the first element has index 0.

For a one-dimensional array:

import numpy as np

a = np.array([10, 20, 30, 40, 50])
print(a[2])     # 30
print(a[1:4])   # [20 30 40]
print(a[::-1])  # reversed array

For a two-dimensional array, indices are specified as [row, column]:

m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
print(m[1, 2])      # 6
print(m[0, :])      # first row
print(m[:, 1])      # second column
print(m[0:2, 1:3])  # selected subarray

The general slicing form is $start:stop:step$, where stop is excluded. NumPy slicing commonly returns a view rather than an independent copy, so changes to a slice may affect the original array.

The important attributes of a NumPy array describe its size, structure, and data representation.

For an array a, the major attributes are:

• a.ndim: number of dimensions or axes.
• a.shape: size of the array along each axis.
• a.size: total number of elements.
• a.dtype: data type of the elements.
• a.itemsize: number of bytes occupied by one element.
• a.nbytes: total memory used by the elements.

Example:

import numpy as np

a = np.array([[1, 2, 3], [4, 5, 6]], dtype=np.int32)
print(a.ndim)     # 2
print(a.shape)    # (2, 3)
print(a.size)     # 6
print(a.dtype)    # int32
print(a.itemsize) # 4
print(a.nbytes)   # 24

These attributes are useful for checking whether arrays have compatible shapes and data types before performing operations.

NumPy performs arithmetic operations element by element when arrays have compatible shapes.

For arrays and of the same shape:

The same principle applies to subtraction, multiplication, division, and exponentiation.

import numpy as np

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

print(A + B)
print(A * B)  # element-wise multiplication
print(A ** 2) # element-wise square

Element-wise multiplication is different from matrix multiplication. Matrix multiplication is defined by:

In NumPy, matrix multiplication can be performed using A @ B or np.matmul(A, B). Thus, A * B and A @ B generally produce different results.

Broadcasting allows NumPy to perform operations on arrays with different but compatible shapes without explicitly replicating data.

NumPy compares dimensions from right to left. Two dimensions are compatible when:

• They are equal, or
• One of them is equal to 1.

Missing dimensions are treated as having size 1.

Example:

import numpy as np

A = np.array([[1, 2, 3], [4, 5, 6]])
b = np.array([10, 20, 30])
print(A + b)

The shape of A is (2, 3) and the shape of b is (3,). NumPy conceptually applies b to each row, producing:

Broadcasting improves performance and avoids unnecessary copies, but incompatible shapes cause a broadcasting error.

NumPy provides several functions for changing the structure of arrays.

• reshape() changes the shape without changing the data. The total number of elements must remain the same.
• ravel() returns a flattened one-dimensional view whenever possible.
• flatten() returns a flattened copy of the array.
• resize() changes the shape and may add or remove elements.

Example:

import numpy as np

a = np.array([[1, 2, 3], [4, 5, 6]])
print(a.reshape(3, 2))
print(a.ravel())
print(a.flatten())

b = np.array([1, 2, 3, 4])
b.resize(2, 3)
print(b)

The distinction between a view and a copy is important. Changes made through a view may affect the original array, whereas changes to a copy do not. The product of the new dimensions in reshape() must equal the original number of elements.

NumPy array manipulation functions change the arrangement or combination of arrays.

• transpose() exchanges axes. For a two-dimensional matrix, it changes rows into columns.
• concatenate() joins arrays along an existing axis.
• vstack() joins arrays vertically, by rows.
• hstack() joins arrays horizontally, by columns.
• split() divides an array into multiple subarrays.

Example:

import numpy as np

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

print(A.T)
print(np.vstack((A, B)))
print(np.hstack((A, np.array([[7], [8]]))))
print(np.split(np.array([1, 2, 3, 4]), 2))

When using concatenate(), the dimensions along all non-concatenated axes must match. These functions are useful for constructing larger matrices and separating data into manageable sections.

Reduction functions combine multiple array elements into summary values. Common functions include:

• np.sum() computes the total.
• np.mean() computes the arithmetic average.
• np.min() and np.max() find the smallest and largest values.
• np.std() calculates standard deviation.
• np.argmin() and np.argmax() return positions of extreme values.

For a matrix, the axis argument specifies the direction of reduction:

import numpy as np

A = np.array([[1, 2, 3], [4, 5, 6]])
print(A.sum())       # total of all elements
print(A.sum(axis=0)) # column-wise sums
print(A.sum(axis=1)) # row-wise sums
print(A.mean(axis=0))

For a two-dimensional array, axis=0 reduces rows and produces one result per column, while axis=1 reduces columns and produces one result per row. The keepdims=True option preserves the number of dimensions.

For the matrix

its determinant is:

If , the matrix is nonsingular and its inverse is:

In NumPy:

import numpy as np

A = np.array([[2, 1], [1, 3]], dtype=float)
determinant = np.linalg.det(A)
inverse = np.linalg.inv(A)
print(determinant)
print(inverse)

The inverse satisfies , where is the identity matrix. In numerical applications, solving a system with np.linalg.solve(A, b) is generally preferred to explicitly calculating the inverse.

Matrix decomposition expresses a matrix as a product of simpler matrices. It is useful for solving systems, calculating least-squares solutions, and analyzing numerical data.

• LU decomposition:
  
  where is lower triangular and is upper triangular. It is useful for solving several systems with the same coefficient matrix.

• QR decomposition:
  
  where has orthonormal columns and is upper triangular. It is commonly used for least-squares problems.

• Singular value decomposition:
  
  where and are orthogonal matrices and contains singular values. SVD is useful in dimensionality reduction, data compression, and determining matrix rank.

SciPy provides routines such as scipy.linalg.lu(), scipy.linalg.qr(), and scipy.linalg.svd(). Each decomposition has different numerical and computational advantages.

A system of linear equations can be written as:

where is the coefficient matrix, is the vector of unknowns, and is the constant vector. If is nonsingular, the system has the unique solution , although explicitly computing the inverse is usually inefficient.

NumPy provides np.linalg.solve():

import numpy as np

A = np.array([[2, 1], [1, 3]], dtype=float)
b = np.array([5, 6], dtype=float)
x = np.linalg.solve(A, b)
print(x)

The result can be verified using:

print(np.allclose(A @ x, b))

solve() is generally more accurate and efficient than calculating np.linalg.inv(A) @ b. If the matrix is singular or not square, other techniques such as least-squares methods may be required.

A sparse matrix contains mostly zero elements, whereas a dense matrix contains a relatively large proportion of nonzero elements.

For example, a matrix of size has positions. If only a few thousand positions are nonzero, storing every zero is wasteful.

Sparse matrices are useful because they:

• Store only nonzero values and their locations.
• Require less memory.
• Reduce computation for operations that involve mostly zeros.
• Are common in graphs, networks, text data, finite-element models, and recommendation systems.

A dense representation stores every element in a rectangular block. A sparse representation stores values together with structural information such as row indices, column indices, or row pointers. The best sparse format depends on the operations to be performed.

SciPy supports several sparse matrix formats.

• COO (Coordinate format): stores three arrays: data, row, and col. Each nonzero value is represented by its value and its row and column indices. COO is convenient for constructing a sparse matrix.
• CSR (Compressed Sparse Row): stores data, indices, and indptr. The indptr array identifies where each row starts and ends. CSR is efficient for row slicing and matrix-vector multiplication.
• CSC (Compressed Sparse Column): is the column-oriented counterpart of CSR. It is efficient for column slicing and column-based operations.

Example:

from scipy.sparse import coo_matrix

values = [5, 8, 3]
rows = [0, 1, 2]
cols = [1, 2, 0]
S = coo_matrix((values, (rows, cols)), shape=(3, 3))
print(S.toarray())
print(S.tocsr())

COO is mainly a construction format, while CSR and CSC are generally better for arithmetic and access operations.

SciPy is an open-source Python library for scientific and technical computing. It is built on NumPy and provides specialized modules for linear algebra, optimization, integration, statistics, signal processing, and sparse computation.

The scipy.linalg module provides functions for:

• Matrix inverse and determinant.
• Eigenvalues and eigenvectors.
• LU, QR, and Schur decompositions.
• Matrix exponentials and logarithms.
• Solving linear systems.

Example:

import numpy as np
from scipy import linalg

A = np.array([[2, 1], [1, 2]], dtype=float)
w, v = linalg.eig(A)
print(w)
print(v)

Although NumPy contains many basic linear-algebra functions, SciPy often provides a broader and more specialized collection of routines, with additional numerical algorithms and sparse alternatives.

For a square matrix , a nonzero vector is an eigenvector if:

where is the corresponding eigenvalue. Rearranging gives the characteristic equation:

The roots of this equation are eigenvalues, and the corresponding vectors are eigenvectors.

Example using SciPy:

import numpy as np
from scipy.linalg import eig

A = np.array([[4, 1], [2, 3]], dtype=float)
eigenvalues, eigenvectors = eig(A)
print(eigenvalues)
print(eigenvectors)

Each column of eigenvectors corresponds to an eigenvalue. Eigenvalues and eigenvectors are used in stability analysis, principal component analysis, vibration analysis, and differential equations.

Dense matrices store all entries, including zeros, while sparse matrices store only nonzero entries and structural information.

A sparse representation should generally be preferred when:

• The matrix has a high proportion of zeros.
• The matrix is too large for practical dense storage.
• Operations such as matrix-vector multiplication can exploit the zero pattern.
• The application naturally produces sparse data, such as graph adjacency matrices.

Example:

import numpy as np
from scipy.sparse import csr_matrix

A = csr_matrix([[0, 5, 0], [0, 0, 2], [7, 0, 0]])
x = np.array([1, 2, 3])
print(A @ x)

For a small or nearly full matrix, dense storage may be faster and simpler. Converting a very large sparse matrix to a dense matrix can consume excessive memory, so operations should usually remain in sparse form.

Vectorization means expressing an operation on an entire array rather than processing elements one at a time in a Python loop.

Non-vectorized code:

result = []
for value in values:
    result.append(value ** 2 + 1)

Vectorized code:

import numpy as np

values = np.array([1, 2, 3, 4])
result = values ** 2 + 1

Advantages include:

• Faster execution because operations are performed in optimized compiled code.
• Shorter and clearer programs.
• Better use of low-level numerical libraries and hardware.
• Easy combination with broadcasting and universal functions.

Vectorization does not mean that loops are absent internally; rather, the loops are moved from slow Python-level execution into optimized library implementations.

A universal function, or ufunc, is a NumPy function that operates element by element on arrays. Ufuncs support vectorization, broadcasting, and usually efficient compiled implementations.

Examples include:

• np.sqrt() for square roots.
• np.exp() for exponentials.
• np.sin() and np.cos() for trigonometric functions.
• np.add(), np.multiply(), and np.power() for arithmetic.

Example:

import numpy as np

x = np.array([0, np.pi / 2, np.pi])
y = np.sin(x)
print(y)

The expression np.sin(x) applies the sine function to every element without an explicit loop. Ufuncs can also accept scalar values, arrays, and broadcast-compatible arrays. They are an important foundation of fast vectorized numerical programming in Python.

Suppose two arrays contain corresponding points:

The Euclidean distance between corresponding points is:

A vectorized implementation is:

import numpy as np

P = np.array([[1, 2], [3, 4], [5, 6]], dtype=float)
Q = np.array([[2, 4], [1, 1], [6, 8]], dtype=float)

distances = np.sqrt(np.sum((P - Q) ** 2, axis=1))
print(distances)

Explanation:

• P - Q computes coordinate differences element by element.
• ** 2 squares every difference.
• np.sum(..., axis=1) adds the squared coordinate differences for each point.
• np.sqrt() computes the square root of each sum.

This approach is concise and faster than a Python loop for large arrays.