NumPy Training

NumPy is a foundational Python library for numerical computing, providing fast, efficient support for multi‑dimensional arrays and mathematical operations. It originated in the 1990s from earlier projects (Numeric and Numarray) and was unified into NumPy in 2005 by Travis Oliphant. It’s useful whenever you need high‑performance math, data analysis, or scientific computing, and its features include array objects, broadcasting, linear algebra, random number generation, and interoperability with other libraries.

What is NumPy?

• Definition: NumPy (short for Numerical Python) is an open‑source library that adds support for large, multi‑dimensional arrays and matrices, along with a wide collection of mathematical functions.
• Core Idea: Replace slow Python loops with fast, vectorized operations written in C under the hood.
• Role: It’s the backbone of the Python scientific ecosystem — libraries like Pandas, SciPy, scikit‑learn, TensorFlow, and PyTorch all rely on NumPy arrays.

📜 History

• 1995: Numeric created by Jim Hugunin as the first array package for Python.
• Early 2000s: Numarray developed as a more flexible replacement, but both Numeric and Numarray had trade‑offs.
• 2005: Travis Oliphant unified the two projects into NumPy, adding extensive modifications.
• Community Growth: NumPy became the standard for numerical computing in Python, supported by NumFOCUS and widely adopted in science and engineering.
• Impact: Used in major scientific breakthroughs like gravitational wave analysis and black hole imaging

Why and When is NumPy Useful?

When: Anytime you need fast numerical operations on large datasets, especially in data science, machine learning, physics, engineering, or finance.

Why:

• Pure Python is slow for math on large lists.
• NumPy arrays are stored in contiguous memory and processed with optimized C routines.
• Enables vectorization (operations on entire arrays at once).
• Provides interoperability with other scientific libraries.

Key Features

✅ Key Takeaway

NumPy is the engine of numerical computing in Python. Born from early array packages in the 1990s, it became the unified standard in 2005 and now underpins nearly all modern data science and AI libraries. It’s most useful when working with large datasets, mathematical modeling, or scientific computing, thanks to its speed, broadcasting, and rich set of features.

Array Types in NumPy

In NumPy, arrays are homogeneous — all elements share the same data type (dtype). NumPy supports a wide range of numerical, string, and object types. You can specify the type when creating an array, check it with .dtype, and convert between types using .astype().

NumPy supports many more types than standard Python. The main categories are:

You can also use character codes:

• 'i' → integer, 'u' → unsigned integer, 'f' → float, 'c' → complex, 'S' → string, 'U' → unicode

Creating Arrays with Specific Types

import numpy as np

# Integer array
arr1 = np.array([1, 2, 3], dtype=np.int32)

# Float array
arr2 = np.array([1, 2, 3], dtype=np.float64)

# String array
arr3 = np.array([1, 2, 3], dtype='U')  # Unicode strings

Check type:

print(arr1.dtype)  # int32
print(arr2.dtype)  # float64
print(arr3.dtype)  # <U11 (unicode string of length 11)

Conversions Between Types

The main method is .astype():

arr = np.array([10, 20, 30], dtype=np.int32)

# Convert to float
arr_float = arr.astype(np.float64)

# Convert to complex
arr_complex = arr.astype(np.complex128)

print(arr_float.dtype)   # float64
print(arr_complex.dtype) # complex128

• .astype() creates a new array with the desired type.
• If conversion is impossible (e.g., converting 'apple' to integer), NumPy raises a ValueError

Multidimensional arrays

Multidimensional arrays in NumPy (called ndarray) are powerful data structures that allow you to store and manipulate data in multiple dimensions efficiently. They are widely used for mathematical computations, image processing, and scientific analysis because they support fast operations, slicing, reshaping, and broadcasting. Multidimensional arrays are a core feature of NumPy (Numerical Python), which is the fundamental package for scientific computing in Python. They are primarily handled by the ndarray object.

ndarray: The Core Object

The NumPy ndarray (N-dimensional array) is a grid of values, all of the same type, and is indexed by a tuple of non-negative integers.

• 1D Array (Vector): A single row or column of elements.
• 2D Array (Matrix): A grid with rows and columns.
• 3D Array: An array of matrices (like a cube or layers).
• nD Array: Generalization to any number of dimensions.

Every NumPy array has several important attributes that describe its structure:

Creating Multidimensional Arrays

You can create these arrays in several ways:

From Python Lists/Tuples:

import numpy as np

# 1D Array
a = np.array([1, 2, 3])
# a.shape is (3,)

# 2D Array (Matrix)
b = np.array([[1, 2, 3], [4, 5, 6]])
# b.shape is (2, 3) -> 2 rows, 3 columns

Using built-in NumPy functions:

# An array of zeros with a specific shape
c = np.zeros((2, 3, 4))
# c.shape is (2, 3, 4) -> 2 blocks, 3 rows, 4 columns

# An array of ones
d = np.ones((5, 5))

# A range of numbers reshaped into a 2D array
e = np.arange(12).reshape((3, 4))
# e will be a 3x4 matrix with values 0 to 11

print(e[2, 3])

Creating arrays from lists and other Python structures

NumPy makes it very easy to create arrays from lists, tuples, ranges, and other Python structures. This is often the first step when moving from plain Python into numerical computing.

import numpy as np

# From a simple list
arr1 = np.array([1, 2, 3, 4])
print(arr1)        # [1 2 3 4]
print(arr1.shape)  # (4,)

# From Nested Lists (Multidimensional)
arr2 = np.array([[1, 2, 3], [4, 5, 6]])
print(arr2)
# [[1 2 3]
#  [4 5 6]]
print(arr2.shape)  # (2, 3)

# From Tuples
arr_tuple = np.array((7, 8, 9))
print(arr_tuple)  # [7 8 9]

# From Ranges
arr_range = np.array(range(5))
print(arr_range)  # [0 1 2 3 4]

# From Comprehensions or Generators
arr_gen = np.array([x**2 for x in range(5)])
print(arr_gen)  # [0 1 4 9 16]

• Homogeneity: All elements must be of the same type. If mixed, NumPy will upcast (e.g., integers + floats → floats).
• Efficiency: Arrays are stored in contiguous memory blocks, making them faster than lists.
• Shape inference: NumPy figures out dimensions automatically based on nesting depth.

Common Intrinsic Array Creation Functions

1. np.zeros()

Creates an array filled with zeros.

import numpy as np

arr = np.zeros((3, 4))
print(arr)
# [[0. 0. 0. 0.]
#  [0. 0. 0. 0.]
#  [0. 0. 0. 0.]]

• Shape (3,4) → 3 rows × 4 columns.
• Default dtype is float.

2. np.ones()

Creates an array filled with ones.

arr = np.ones((2, 3), dtype=int)
print(arr)
# [[1 1 1]
#  [1 1 1]]

3. np.empty()

Creates an array without initializing values (fast, but contents are random garbage).

arr = np.empty((2, 2))
print(arr)
# [[1.23e-322 0.00e+000]
#  [ ... random values ... ]]

4. np.arange()

Creates evenly spaced values within a range (like Python’s range, but returns an array).

arr = np.arange(0, 10, 2)
print(arr)  # [0 2 4 6 8]

Parameters: start, stop, step.

5. np.linspace()

Creates evenly spaced numbers over a specified interval.

arr = np.linspace(0, 1, 5)
print(arr)  # [0.   0.25 0.5  0.75 1. ]

Parameters: start, stop, num (number of points).

6. np.eye()

Creates an identity matrix.

arr = np.eye(3)
print(arr)
# [[1. 0. 0.]
#  [0. 1. 0.]
#  [0. 0. 1.]]

7. np.full()

Creates an array filled with a specific value.

arr = np.full((2, 3), 7)
print(arr)
# [[7 7 7]
#  [7 7 7]]

add, remove, and sort elements in NumPy arrays

Adding Elements

NumPy arrays don’t support direct append like lists, but you can use np.append or np.insert.

import numpy as np

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

# Append at the end
arr_add = np.append(arr, 4)
print(arr_add)   # [1 2 3 4]

# Insert at index 1
arr_insert = np.insert(arr, 1, 99)
print(arr_insert)  # [ 1 99  2  3]

• np.append(arr, values) → adds values at the end.
• np.insert(arr, index, values) → inserts values at a specific position.
• Note: These return new arrays, since NumPy arrays are immutable in size.

Removing Elements

Use np.delete to remove elements by index.

arr = np.array([10, 20, 30, 40])

# Remove element at index 2
arr_del = np.delete(arr, 2)
print(arr_del)   # [10 20 40]

# Remove multiple indices
arr_del_multi = np.delete(arr, [0, 3])
print(arr_del_multi)  # [20 30]

Sorting Elements

NumPy provides np.sort and ndarray.sort.

arr = np.array([3, 1, 4, 2])

# Sorted copy
sorted_arr = np.sort(arr)
print(sorted_arr)  # [1 2 3 4]

# In-place sort
arr.sort()
print(arr)  # [1 2 3 4]

Sorting along axes (multidimensional arrays)

arr2d = np.array([[3, 1, 2],
                  [9, 7, 8]])

# Sort each row
print(np.sort(arr2d, axis=1))
# [[1 2 3]
#  [7 8 9]]

# Sort each column
print(np.sort(arr2d, axis=0))
# [[3 1 2]
#  [9 7 8]]

copies vs views

What’s the difference?

• Copy: A completely new array object with its own data. Changes to the copy do not affect the original.
• View: A new array object that looks at the same data buffer as the original. Changes to the view do affect the original (and vice versa).

1. Copy

import numpy as np

arr = np.array([1, 2, 3, 4])
copy_arr = arr.copy()

copy_arr[0] = 99

print("Original:", arr)   # [1 2 3 4]
print("Copy:", copy_arr)  # [99 2 3 4]

• arr.copy() creates a deep copy.
• Modifying copy_arr does not change arr.

2. View (slicing)

arr = np.array([1, 2, 3, 4])
view_arr = arr[1:3]   # slice creates a view

view_arr[0] = 99

print("Original:", arr)   # [1 99 3 4]
print("View:", view_arr)  # [99 3]

• Slicing returns a view.
• Changing view_arr also changes arr.

3. View (reshape)

arr = np.array([1, 2, 3, 4, 5, 6])
reshaped = arr.reshape(2, 3)  # reshape often returns a view

reshaped[0, 0] = 99

print("Original:", arr)     # [99 2 3 4 5 6]
print("Reshaped:", reshaped)
# [[99  2  3]
#  [ 4  5  6]]

• reshape usually returns a view (if possible).
• Both arrays share the same data buffer.

4. View

• A view is a new array object that looks at the same data as the original.
• Created automatically by operations like slicing or reshaping (when possible).
• Changes in the view affect the original array.

import numpy as np

arr = np.array([1, 2, 3, 4])
view_arr = arr.view()

view_arr[1] = 5
print('arr : ', arr) # [1 5 3 4]
print('view_arr : ', view_arr) # [1 5 3 4]

5. ravel()

• Returns a flattened view of the array (1D).
• If possible, it avoids copying data (so changes may affect the original).
• Faster and memory‑efficient compared to flatten().

arr = np.array([[1, 2], [3, 4]])
r = arr.ravel()
# or r= np.reshape(arr, -1)
print(r)   # [1 2 3 4]

r[0] = 99
print(arr) # [[99  2]
           #  [ 3  4]]

6. flatten()

• Returns a flattened copy of the array (always 1D).
• Independent of the original — changes won’t affect the source.

arr = np.array([[1, 2], [3, 4]])
f = arr.flatten()

print(f)   # [1 2 3 4]

f[0] = 99
print(arr) # [[1 2]
           #  [3 4]]   # unchanged

What is .base?

• Every NumPy array has a .base attribute.
• It tells you whether the array owns its data or is a view of another array.
• If .base is None, the array owns its data (it’s independent).
• If .base is another array, it means the current array is a view into that other array’s data buffer.

1. Copy (independent)

import numpy as np

arr = np.array([1, 2, 3, 4])
copy_arr = arr.copy()

print(copy_arr.base)   # None → owns its data

2. View (slice)

arr = np.array([1, 2, 3, 4])
view_arr = arr[1:3]

print(view_arr.base is arr)  # True → shares data with arr

3. Reshape (view if possible)

arr = np.arange(6)
reshaped = arr.reshape(2, 3)

print(reshaped.base is arr)  # True → reshaped is a view of arr

Why it matters

• .base helps you check if two arrays share memory.
• Useful for debugging unexpected changes when modifying slices or reshaped arrays.
• Helps you decide whether you need a copy to avoid side effects.

Rule of thumb:

• If .base is None, the array is independent.
• If .base points to another array, it’s a view and shares data.

Indexing in NumPy

Basic Indexing

Single element indexing works like Python lists:

import numpy as np
arr = np.array([10, 20, 30, 40])
print(arr[0])   # 10
print(arr[-1])  # 40

Multidimensional arrays use commas to separate dimensions:

arr2d = np.array([[1, 2, 3],
                  [4, 5, 6]])
print(arr2d[0, 1])  # 2
print(arr2d[1, 2])  # 6

Slicing

Syntax: [start:stop:step]

Includes start, excludes stop.

Defaults: start=0, stop=length, step=1.

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

2D slicing

arr2d = np.array([[1, 2, 3],
                  [4, 5, 6],
                  [7, 8, 9]])
print(arr2d[:2, 1:])  
# [[2 3]
#  [5 6]]

Advanced Indexing

Fancy Indexing

Use lists/arrays of indices:

arr = np.array([10, 20, 30, 40, 50])
print(arr[[0, 2, 4]])  # [10 30 50]

2D fancy indexing

arr2d = np.array([[1, 2], [3, 4], [5, 6]])
print(arr2d[[0, 2], [1, 0]])   # [2 5]

Boolean Indexing

Use conditions to filter:

arr = np.array([10, 20, 30, 40, 50])
print(arr[arr > 25])   # [30 40 50]

• Slicing returns a view → changes affect the original array.
• Fancy/boolean indexing returns a copy → changes don’t affect the original.
• Indexing and slicing are optimized in NumPy, making them much faster than Python lists for large datasets

join or split

In NumPy, you often need to join (concatenate/stack) arrays together or split (partition) them into smaller arrays. Let’s go through the main tools:

Joining Arrays

1. np.concatenate

Joins arrays along an existing axis.

The most general one

• Joins arrays along an existing axis
• Does not increase dimensions

import numpy as np

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

print(np.concatenate((a, b)))   # [1 2 3 4 5 6]

arr2d = np.array([[3, 1, 7],
                  [8, 9, 2]])
arr2d1 = np.array([[3, 2, 7],
                  [1, 9, 2]])
print(np.concatenate((arr2d, arr2d1), axis=0)) #or np.concatenate((arr2d, arr2d1))
# array([[3, 1, 7],
#       [8, 9, 2],
#       [3, 2, 7],
#       [1, 9, 2]])

print(np.concatenate((arr2d, arr2d1), axis=1))
#array([[3, 1, 7, 3, 2, 7],
#      [8, 9, 2, 1, 9, 2]])

2. np.stack

Joins arrays along a new axis.

The only one that creates a new axis

• All arrays must have the same shape
• Output has one more dimension

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

print(np.stack((a, b), axis=0))
# [[1 2 3]
#  [4 5 6]]

print(np.stack((a, b), axis=1))
# [[1 4]
# [2 5]
# [3 6]]

arr2d = np.array([[3, 1, 7],
                  [8, 9, 2]])
arr2d1 = np.array([[3, 2, 7],
                  [1, 9, 2]])


print(np.stack((arr2d, arr2d1), axis=0))
#array([[[3, 1, 7],
#       [8, 9, 2]],

#      [[3, 2, 7],
#       [1, 9, 2]]])

print(np.stack((arr2d, arr2d1), axis=1))
# array([[[3, 1, 7],
#       [3, 2, 7]],

#      [[8, 9, 2],
#       [1, 9, 2]]])

3. np.hstack and np.vstack

Horizontal and vertical stacking.

Vertical stack (row-wise)

Basically:

np.vstack([a, b]) == np.concatenate([a, b], axis=0)

For 1D arrays, vstack first turns them into rows.

Horizontal stack (column-wise)

• For 2D+ arrays:

np.hstack([a, b]) == np.concatenate([a, b], axis=1)

Important gotcha for 1D arrays. For 1D arrays, hstack behaves like concatenate(axis=0).

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

print(np.hstack((a, b)))   # [1 2 3 4 5 6]
print(np.vstack((a, b)))
# [[1 2 3]
#  [4 5 6]]

arr2d = np.array([[3, 1, 7],
                  [8, 9, 2]])
arr2d1 = np.array([[3, 2, 7],
                  [1, 9, 2]])

print(np.hstack((arr2d, arr2d1)))
#array([[3, 1, 7, 3, 2, 7],
#       [8, 9, 2, 1, 9, 2]])

print(np.vstack((arr2d, arr2d1)))
#array([[3, 1, 7],
#       [8, 9, 2],
#       [3, 2, 7],
#       [1, 9, 2]])

• Need a new dimension? → np.stack
• Just joining along an axis? → np.concatenate
• Stack rows? → np.vstack
• Stack columns? → np.hstack

4. np.column_stack

Stacks 1D arrays as columns into a 2D array.

Stack arrays as columns

• For 1D arrays → turns them into columns first
• For 2D arrays → behaves like hstack
• Does NOT create a new axis

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

print(np.column_stack((a, b)))
# [[1 4]
#  [2 5]
#  [3 6]]

arr2d = np.array([[3, 1, 7],
                  [8, 9, 2]])
arr2d1 = np.array([[3, 2, 7],
                  [1, 9, 2]])

print(np.column_stack([arr2d, arr2d1]))
# array([[3, 1, 7, 3, 2, 7],
#       [8, 9, 2, 1, 9, 2]])

Splitting Arrays

NumPy has a family of split functions that all do the same basic job — cut an array into smaller arrays — but they differ in how strict they are and which axis they operate on by default.

1. np.split

Splits an array into multiple and equal-sized subarrays.

Key rule

The split must be even, or it raises an error.

np.split(arr, indices_or_sections, axis=0)

When to use

• You know the array divides evenly
• You want strict correctness

import numpy as np

arr = np.array([1, 2, 3, 4, 5, 6])
print(np.split(arr, 3))
# [array([1, 2]), array([3, 4]), array([5, 6])]

np.split(arr, 4)
# ValueError: array split does not result in an equal division

A = np.array([[ 1,  2,  3,  4],
              [ 5,  6,  7,  8],
              [ 9, 10, 11, 12]])

np.split(A, 3)
#[
# array([[1, 2, 3, 4]]),
# array([[5, 6, 7, 8]]),
# array([[ 9, 10, 11, 12]])
# ]

np.split(A, 2, axis=1)

#[
# array([[1, 2],
#        [5, 6],
#        [9, 10]]),
# array([[3, 4],
#        [7, 8],
#        [11, 12]])
#]

np.split(A, 2, axis=0)
# ValueError: array split does not result in an equal division

2. np.array_split

Like split, but allows unequal splits. Splits an array into approximately equal-sized chunks.

Key difference from np.split

Works even when the array does not divide evenly.

arr = np.array([1, 2, 3, 4, 5])
print(np.array_split(arr, 2))
# [array([1, 2, 3]), array([4, 5])]

np.array_split(A, 2, axis=0)
#[
# array([[1, 2, 3, 4],
#        [5, 6, 7, 8]]),
# array([[ 9, 10, 11, 12]])
#]


np.array_split(A, 3, axis=1)
#[
# array([[1, 2],
#        [5, 6],
#        [9, 10]]),
# array([[3],
#        [7],
#        [11]]),
# array([[4],
#        [8],
#        [12]])
#]

How it handles leftovers

• Earlier chunks get the extra elements

When to use

• You don’t know the exact size ahead of time
• You want a “best effort” split

3. np.hsplit and np.vsplit

np.hsplit

Splits an array horizontally (along columns) Equivalent to

np.split(arr, ..., axis=1)

When to use

• You’re thinking in terms of columns
• You want readability (hsplit is very explicit)

np.vsplit

Splits an array vertically (along rows). Equivalent to

np.split(arr, ..., axis=0)

Works on

• 2D arrays (or higher)

When to use

• You’re thinking in terms of rows
• Clearer than specifying axis=0

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

print(np.hsplit(arr2d, 3))  # Split into 3 columns
# [array([[1],[4]]), array([[2],[5]]), array([[3],[6]])]

print(np.vsplit(arr2d, 2))  # Split into 2 rows
# [array([[1,2,3]]), array([[4,5,6]])]

np.hsplit(A, 2)  # np.split(A, 2, axis=1)
#[
# array([[1, 2],
#        [5, 6],
#        [9, 10]]),
# array([[3, 4],
#        [7, 8],
#        [11, 12]])
#]

np.hsplit(A, 3) # Error (same rule as split).

np.vsplit(A, 3)  # np.split(A, 3, axis=0)
#[
# array([[1, 2, 3, 4]]),
# array([[5, 6, 7, 8]]),
# array([[ 9, 10, 11, 12]])
#]

np.vsplit(A, 2) # Error (rows don’t divide evenly).

Mental model

• split → “I demand equality”
• array_split → “Do your best”
• hsplit → “Cut left to right”
• vsplit → “Cut top to bottom”

Arithmetic operations and functions

NumPy is designed for fast, vectorized arithmetic. You can perform operations element‑wise across arrays, and it also provides a rich set of mathematical functions. Let’s break it down:

Basic Arithmetic Operations

All operations are element‑wise (applied to each element of the array).

import numpy as np

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

print(a + b)   # [5 7 9]   (addition)
print(a - b)   # [-3 -3 -3] (subtraction)
print(a * b)   # [4 10 18] (multiplication)
print(a / b)   # [0.25 0.4 0.5] (division)
print(a ** 2)  # [1 4 9]   (power)

Works for scalars too:

print(a + 10)   # [11 12 13]
print(a * 2)    # [2 4 6]

Arithmetic Functions

NumPy provides explicit functions for basic arithmetic:

import numpy as np

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

print(np.add(a, b))       # [5 7 9]
print(np.subtract(a, b))  # [-3 -3 -3]
print(np.multiply(a, b))  # [4 10 18]
print(np.divide(a, b))    # [0.25 0.4 0.5]
print(np.power(a, 2))     # [1 4 9]
print(np.mod(a, b))       # [1 2 3] (remainder)

These are equivalent to operators:

• np.add(a, b) ↔ a + b
• np.multiply(a, b) ↔ a * b
• np.divide(a, b) ↔ a / b

Matrix Operations

NumPy distinguishes between element‑wise multiplication and matrix multiplication.

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

print(A * B)        # element-wise → [[5 12] [21 32]]
print(A @ B)        # matrix multiplication → [[19 22] [43 50]]
print(np.dot(A, B)) # same as @

Aggregate Functions

Operate across entire arrays or along axes.

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

print(np.sum(arr))        # 21
print(np.sum(arr, axis=0)) # [5 7 9] (column sums)
print(np.sum(arr, axis=1)) # [6 15]  (row sums)
print(np.prod(arr))       # 720
print(np.average(arr))    # 3.5
print(np.std(arr))        # 1.707825127659933
print(np.mean(arr))       # 3.5
print(np.min(arr))        # 1
print(np.max(arr))        # 6

Trigonometric & Exponential Functions

NumPy provides vectorized versions of common math functions.

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

print(np.sin(x))   # [0. 1. 0.]
print(np.cos(x))   # [1. 0. -1.]
print(np.exp(x))   # [1. 4.81 23.14]
print(np.log([1, np.e, np.e**2]))  # [0. 1. 2.]

Comparison & Logical Operations

Element‑wise comparisons return boolean arrays.

a = np.array([1, 2, 3])
b = np.array([2, 2, 2])

print(a > b)   # [False False  True]
print(a == b)  # [False  True False]

# Boolean logic
mask = (a > 1) & (a < 3)
print(mask)    # [False  True False]

Comparison Functions

print(np.equal(a, b))      # [False False False]
print(np.not_equal(a, b))  # [ True  True  True]
print(np.greater(a, b))    # [False False False]
print(np.less(a, b))       # [ True  True  True]

Broadcasting

broadcasting is one of NumPy’s most powerful features. It explains how NumPy handles arithmetic operations between arrays of different shapes without explicitly copying data.

Broadcasting Rules

When operating on two arrays:

1. Compare their shapes from right to left (last dimension first).
2. Dimensions are compatible if:

• They are equal, OR
• One of them is 1.

3. If compatible, NumPy “broadcasts” the smaller array across the larger one.

If not compatible → ValueError.

1. Scalar with array

import numpy as np
arr = np.array([1, 2, 3])
print(arr + 10)   # [11 12 13]

Scalar 10 is broadcast to match shape (3,).

2. 1D with 2D

arr2d = np.array([[1, 2, 3],
                  [4, 5, 6]])
vec = np.array([10, 20, 30])

print(arr2d + vec)

# [[11 22 33]
# [14 25 36]]

vec shape (3,) is broadcast across rows of arr2d shape (2,3).

3. Column vector with row vector

col = np.array([[1],
                [2],
                [3]])   # shape (3,1)
row = np.array([10, 20, 30])  # shape (3,)

print(col + row)

# [[11 21 31] 
# [12 22 32]
# [13 23 33]]

Shapes (3,1) and (3,) → broadcast to (3,3).

4. Incompatible shapes

a = np.array([1, 2, 3])    # shape (3,)
b = np.array([1, 2])       # shape (2,)
print(a + b)   # ❌ ValueError

Shapes (3,) and (2,) are not compatible.

Key Insight

• Broadcasting saves memory and computation by avoiding explicit replication.
• It’s what makes NumPy arithmetic vectorized and efficient.
• Always check shapes with .shape to understand how broadcasting will apply.

How to get unique items and counts

NumPy makes it very easy to find unique elements and their counts in an array. The key function is np.unique.

np.unique

Returns the sorted unique values of an array.

With extra arguments, it can also return:

• Indices of the first occurrence
• Inverse indices (to reconstruct the original array)
• Counts of each unique value

Unique Items Only

import numpy as np

arr = np.array([1, 2, 2, 3, 3, 3, 4])

unique_vals = np.unique(arr)
print(unique_vals)   # [1 2 3 4]

Unique Items with Counts

unique_vals, counts = np.unique(arr, return_counts=True)

print("Unique values:", unique_vals)   # [1 2 3 4]
print("Counts:", counts)               # [1 2 3 1]

Reconstruct Original Array

unique_vals, inverse = np.unique(arr, return_inverse=True)

print("Unique values:", unique_vals)   # [1 2 3 4]
print("Inverse indices:", inverse)     # [0 1 1 2 2 2 3]

# Rebuild original array
print(unique_vals[inverse])            # [1 2 2 3 3 3 4]

Transpose-like operations

1. ndarray.T

The simplest transpose: swaps rows and columns for 2D arrays.

import numpy as np

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

print(arr.T)
# [[1 4]
#  [2 5]
#  [3 6]]

For higher dimensions, .T reverses the order of axes.

2. np.transpose

More general: lets you specify the order of axes.

arr3d = np.arange(24).reshape(2, 3, 4)  # shape (2,3,4)

print(np.transpose(arr3d, (1, 0, 2)).shape)  # (3,2,4)

Here, axes are reordered: axis 0 → 1, axis 1 → 0, axis 2 → 2.

3. np.swapaxes

Swaps two specific axes.

arr3d = np.arange(24).reshape(2, 3, 4)

print(np.swapaxes(arr3d, 0, 1).shape)  # (3,2,4)

Equivalent to np.transpose(arr3d, (1,0,2)).

4. np.moveaxis

Moves one axis to a new position, shifting others accordingly.

arr3d = np.arange(24).reshape(2, 3, 4)

print(np.moveaxis(arr3d, 0, -1).shape)  # (3,4,2)

Axis 0 (size 2) moved to the last position.

Reversing an array

Reversing an array in NumPy is pleasantly simple, and there are multiple clean ways to do it depending on what you want to reverse (the whole array, rows, columns, or specific axes). Let’s walk through the most useful techniques.

1. Reverse a 1D Array

The classic Python slice trick works perfectly in NumPy:

import numpy as np

arr = np.array([1, 2, 3, 4, 5])
print(arr[::-1])
# [5 4 3 2 1]

[::-1] means: start → end, step = –1 (go backwards).

2. Reverse a 2D Array (Different Ways)

A. Reverse rows (flip vertically)

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

print(arr[::-1, :])

# [[7 8 9]
# [4 5 6]
# [1 2 3]]

B. Reverse columns (flip horizontally)

arr[:, ::-1]

# [[3 2 1]
# [6 5 4]
# [9 8 7]]

C. Reverse both rows and columns

arr[::-1, ::-1]

# [[9 8 7]
#  [6 5 4]
#  [3 2 1]]

3. Reverse Along a Specific Axis

NumPy provides a dedicated function:

np.flip

np.flip(arr, axis=0)   # reverse rows
# array([[7, 8, 9],
#       [4, 5, 6],
#       [1, 2, 3]])

np.flip(arr, axis=1)   # reverse columns
# array([[3, 2, 1],
#        [6, 5, 4],
#        [9, 8, 7]])

4. Reverse a Multi‑Dimensional Array

For higher‑dimensional arrays:

np.flip(arr, axis=2) # reverse along axis 2

Advanced Numpy

NumPy universal function (ufunc)

numpy.frompyfunc turns any regular Python function into a NumPy universal function (ufunc) — meaning it can operate element‑wise on NumPy arrays with broadcasting support. This is exactly how the NumPy documentation describes it.

It takes:

• a Python function
• the number of inputs (nin)
• the number of outputs (nout)

and returns a ufunc.

A ufunc is a fast, vectorized function that works on arrays element‑by‑element.

Basic syntax

np.frompyfunc(func, nin, nout)

• func → your Python function
• nin → number of input arguments
• nout → number of outputs

The result is a NumPy ufunc object.

Example: turn a Python function into a ufunc

import numpy as np

def myadd(x, y):
    return x + y

uf = np.frompyfunc(myadd, 2, 1)

print(uf([1, 2, 3], [4, 5, 6]))  # ['5' '7' '9']

Note: outputs are strings because frompyfunc always returns Python objects.

Why use frompyfunc?

• To vectorize any Python function
• To add broadcasting support automatically
• To apply custom logic to NumPy arrays
• To create your own ufuncs when NumPy doesn’t provide one

The resulting ufunc is not fast like built‑in NumPy ufuncs because it still calls Python code. For performance‑critical code, prefer numpy.vectorize or rewrite using NumPy operations. frompyfunc always returns arrays of dtype=object.

strides

Strides in NumPy describe how many bytes you need to jump in memory to move to the next element along each axis of an array. This definition comes directly from NumPy’s documentation, which states that strides are a “tuple of bytes to step in each dimension when traversing an array.”

A NumPy array is stored as a contiguous block of memory. Strides tell NumPy how to interpret that memory:

• Each stride value = number of bytes to skip to move 1 step along that axis.
• The tuple length = number of dimensions.

Example from the documentation: A 2×5 array of 32‑bit integers (4 bytes each) is stored as 40 bytes in memory. Its strides might be:

• First stride: bytes to jump to the next row
• Second stride: bytes to jump to the next column

This matches the official explanation.

import numpy as np

a = np.arange(12).reshape(3, 4)
print(a.strides)

If the dtype is float64 (8 bytes), output might be:

(32, 8)

Meaning:

• Move 32 bytes to go from row i to row i+1
• Move 8 bytes to go from column j to j+1

Why strides matter

Get Marzieh Vahdat’s stories in your inbox

Join Medium for free to get updates from this writer.

Remember me for faster sign in

Strides allow NumPy to:

• Create views without copying data
• Slice arrays efficiently
• Reshape arrays cheaply
• Perform fast vectorized operations

They are essential for understanding memory layout and performance, as emphasized in multiple tutorials.

%timeit

%timeit is an IPython/Jupyter magic command that measures how long a piece of Python code takes to run. It’s one of the most useful tools for quick performance testing.

It doesn’t work in plain Python scripts — only in environments like Jupyter Notebook, IPython shell, or Colab.

It runs your code many times, computes the average execution time, and reports a stable, reliable timing.

You don’t need to choose the number of loops — %timeit automatically picks a good number to get accurate results.

Time a single line:

%timeit x = sum(range(1000))

# 6.03 μs ± 58 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

Time a multi‑line block:

Use the cell magic %%timeit:

%%timeit
x = 0
for i in range(1000):
    x += i

# 16.2 μs ± 72.3 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

Meaning:

1. 16.2 μs per loop

• Your code took 16.2 microseconds on average to run once.
• μs = microseconds = 10−6 seconds.

2. ± 72.3 ns

• The variation (standard deviation) between runs was 72.3 nanoseconds.
• ns = nanoseconds = 10−9 seconds.
• A small deviation means the timing was very stable.

3. 7 runs

• %timeit repeated the whole timing process 7 times.
• It does this to get a reliable average.

4. 100,000 loops each

• In each of the 7 runs, your code was executed 100,000 times.
• %timeit automatically chooses this number to get accurate measurements.

Why %timeit is better than time.time()

• It runs the code many times
• It avoids noise from background processes
• It chooses the number of loops automatically
• It gives statistically meaningful results

Structured arrays

Structured arrays in NumPy let you store tabular, record‑like data inside a single NumPy array, where each element can have multiple named fields with different data types. Think of them as NumPy’s version of a lightweight spreadsheet row or a C struct.

They’re incredibly useful when you want NumPy’s speed but need mixed‑type data.

A structured array has:

• Named fields (like columns)
• Different dtypes per field
• Row‑based records

Example fields:

• "name" → string
• "age" → integer
• "weight" → float

All stored in one array.

Basic example

import numpy as np

dt = np.dtype([('name', 'U10'), ('age', 'i4'), ('weight', 'f4')])

data = np.array([
    ('Alice', 25, 55.5),
    ('Bob',   30, 72.3),
    ('Cara',  22, 48.0)
], dtype=dt)

print(data)

Accessing fields

Get a single column:

data['age']

Get a single row:

data[1]

Get a specific value:

data['weight'][2]

Why structured arrays are useful

• Store mixed types efficiently
• Keep data in a compact binary format
• Fast vectorized operations on individual fields
• Interoperate with binary files, C code, memory‑mapped data
• Great for scientific datasets with many attributes per sample

Defining dtypes

You can define fields with:

np.dtype([('field1', 'i4'), ('field2', 'f8')])

Or nested structures:

np.dtype([('point', [('x', 'f4'), ('y', 'f4')]), ('label', 'U5')])

Example: sorting by a field

np.sort(data, order='age')

A structured array is still a NumPy array, but each element behaves like a record with named attributes.

Nested structured

Nested structured dtypes in NumPy let you build records inside records, very similar to C structs inside structs. A full, realistic example helps show how powerful this can be.

Below is a complete, working example with multiple levels of nesting, field access, and usage.

Imagine modeling a dataset of people, where each person has:

a name (string)

a position with fields:

• x, y, z (floats)

stats with fields:

• height (float)
• weight (float)
• age (int)

This is a perfect case for nested dtypes.

Step 1 — Define nested dtypes

import numpy as np

# Nested dtype for a 3D point
position_dtype = np.dtype([
    ('x', 'f4'),
    ('y', 'f4'),
    ('z', 'f4')
])

# Nested dtype for stats
stats_dtype = np.dtype([
    ('height', 'f4'),
    ('weight', 'f4'),
    ('age', 'i4')
])

# Full dtype combining everything
person_dtype = np.dtype([
    ('name', 'U20'),
    ('position', position_dtype),
    ('stats', stats_dtype)
])

Step 2 — Create an array using the nested dtype

people = np.array([
    ('Alice', (1.0, 2.0, 3.0), (165.0, 55.0, 25)),
    ('Bob',   (4.0, 5.0, 6.0), (180.0, 80.0, 30)),
    ('Cara',  (7.0, 8.0, 9.0), (170.0, 60.0, 22))
], dtype=person_dtype)

Step 3 — Access nested fields

Get all names:

people['name']

Get all x‑coordinates:

people['position']['x']

Get Bob’s weight:

people['stats']['weight'][1]

Get Cara’s full position record:

people['position'][2]

Get Alice’s age:

people['stats']['age'][0]

Step 4 — Modify nested fields

Increase everyone’s z‑coordinate by 10:

people['position']['z'] += 10

Update Bob’s height:

people['stats']['height'][1] = 182.0

Step 5 — Use nested fields in computations

Compute BMI for each person:

height_m = people['stats']['height'] / 100
weight = people['stats']['weight']
bmi = weight / (height_m ** 2)

dates and times

NumPy has a surprisingly rich system for working with dates and times, built around two special dtypes:

• datetime64 → represents absolute points in time
• timedelta64 → represents differences between times

They behave like normal NumPy scalars and arrays, but with powerful vectorized operations.

Let’s walk through the essentials in a clean, practical way.

1. Creating dates (datetime64)

Single date

np.datetime64('2025-01-01')

With time

np.datetime64('2025-01-01T12:30')

With specific units

np.datetime64('2025-01-01', 'D')   # day precision
np.datetime64('2025-01-01T12:30', 'm')  # minute precision

2. Arrays of dates

dates = np.array(['2025-01-01', '2025-01-05', '2025-01-10'], dtype='datetime64[D]')

3. Time differences (timedelta64)

Subtracting two dates gives a timedelta64:

d1 = np.datetime64('2025-01-10')
d2 = np.datetime64('2025-01-01')

delta = d1 - d2
print(delta)        # 9 days
print(delta.dtype)  # timedelta64[D]

4. Adding/subtracting timedeltas

np.datetime64('2025-01-01') + np.timedelta64(10, 'D')
#np.datetime64('2025-01-11')

5. Generating date ranges

NumPy can generate sequences of dates just like arange:

np.arange('2025-01-01', '2025-01-10', dtype='datetime64[D]')

6. Supported time units

NumPy supports many granularities:

• Y → years
• M → months
• W → weeks
• D → days
• h → hours
• m → minutes
• s → seconds
• ms, us, ns, ps, fs, as → sub‑second units

Example:

np.datetime64('2025-01-01T12:00:00', 's')

7. Converting between units

np.datetime64('2025-01-01T12:00') .astype('datetime64[s]')

8. Example: compute ages

birthdays = np.array(['1990-05-10', '2000-01-01'], dtype='datetime64[D]')
today = np.datetime64('2025-01-01')

ages = (today - birthdays) / np.timedelta64(365, 'D')
print(ages)

np.busday_count

np.busday_count is NumPy’s built‑in tool for counting business days (weekdays) between two dates. It’s incredibly handy when you need calendar‑aware calculations without pulling in a full time‑series library like pandas.

Let’s walk through it clearly and practically.

It returns the number of valid business days between two dates:

• Monday–Friday by default
• Excludes weekends
• Can exclude custom holidays
• Can use custom workweek definitions

It does not include the end date.

Basic example

import numpy as np

np.busday_count('2025-01-01', '2025-01-10')

This counts business days from Jan 1 up to (but not including) Jan 10.

Including holidays

holidays = ['2025-01-02', '2025-01-06']

np.busday_count('2025-01-01', '2025-01-10', holidays=holidays)

Any date in holidays is treated as a non‑working day.

Custom workweek

Suppose your workweek is Saturday–Wednesday, and weekends are Thursday–Friday.

np.busday_count(
    '2025-01-01',
    '2025-01-10',
    weekmask='Sat Sun Mon Tue Wed'
)

The weekmask string lists the working days.

Vectorized usage

You can pass arrays:

start = np.array(['2025-01-01', '2025-01-05'], dtype='datetime64[D]')
end   = np.array(['2025-01-10', '2025-01-20'], dtype='datetime64[D]')

np.busday_count(start, end)

NumPy computes each pair element‑wise.

np.is_busday

It returns a Boolean array:

• True → the date is a business day
• False → the date is a weekend or holiday

Decomposition

Eigenvalue Decomposition

Eigenvalue decomposition (or eigendecomposition) is the process of breaking a square matrix into its eigenvalues and eigenvectors, which reveal the matrix’s fundamental behavior. In simple terms: it shows the directions a matrix stretches space and by how much.

For a square matrix A:

Av=λv

• v → an eigenvector (a special direction that doesn’t change orientation under A)
• λ → the eigenvalue (the factor by which A stretches or shrinks v)

If a matrix has enough independent eigenvectors, it can be decomposed as:

• V → matrix of eigenvectors (columns are eigenvectors)
• Λ → diagonal matrix of eigenvalues
• V^−1 → inverse of eigenvector matrix

This is called eigendecomposition.

Step-by-Step Process

1. Find eigenvalues: Solve the characteristic equation

det⁡(A−λI)=0

This polynomial gives possible eigenvalues.

2. Find eigenvectors: For each eigenvalue λ, solve

(A−λI)v=0

to get the corresponding eigenvector(s).

3. Form matrices:

• Put eigenvectors into V (columns).
• Put eigenvalues into diagonal matrix Λ.

4. Check decomposition: Verify A=VΛV^−1.

For square matrices:

Av=λv

vals, vecs = np.linalg.eig(A)

• vals → eigenvalues
• vecs → eigenvectors (columns)

For symmetric matrices, use the more stable version:

vals, vecs = np.linalg.eigh(A)

Example:

import numpy as np

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

vals, vecs = np.linalg.eig(A)

print("Eigenvalues:", vals)
print("Eigenvectors:\n", vecs)

• vals → 1‑D array of eigenvalues
• vecs → 2‑D array whose columns are the eigenvectors

Example: Reconstructing A=VΛV−1

import numpy as np

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

# Eigen decomposition
vals, vecs = np.linalg.eig(A)

# Build Λ (diagonal matrix of eigenvalues)
Lambda = np.diag(vals)

# Reconstruct A
A_reconstructed = vecs @ Lambda @ np.linalg.inv(vecs)

print("Original A:\n", A)
print("Reconstructed A:\n", A_reconstructed)

Applications

• Physics: vibration modes, quantum mechanics
• Machine Learning: PCA (Principal Component Analysis)
• Computer Graphics: transformations, stability analysis
• Finance: risk models, covariance matrices

Singular Value Decomposition (SVD)

SVD is the most general and powerful decomposition. It works for any matrix.

• Works for any matrix (square or rectangular)
• U, V → orthogonal
• Σ → diagonal with singular values

Used in PCA, compression, noise reduction, pseudoinverses, and more.

import numpy as np

A = np.array([[1, 2], [3, 4], [5, 6]])
U, S, Vt = np.linalg.svd(A)

• U → left singular vectors
• S → singular values (1‑D array)
• Vt → right singular vectors (transposed)

import numpy as np

A = np.array([[1, 2], [3, 4], [5, 6]])
U, S, Vt = np.linalg.svd(A)

# Build Sigma with correct shape
Sigma = np.zeros((U.shape[0], Vt.shape[0]))
np.fill_diagonal(Sigma, S)

# Reconstruct A
A_reconstructed = U @ Sigma @ Vt

print("Original A:\n", A)
print("Reconstructed A:\n", A_reconstructed)

Applications of SVD

• Principal Component Analysis (PCA) → dimensionality reduction, data compression
• Latent Semantic Analysis (LSA) → text mining, topic modeling
• Image compression → keep only largest singular values
• Pseudoinverse → solve least‑squares problems
• Noise reduction → discard small singular values

QR Decomposition

QR decomposition (or QR factorization) expresses a matrix as the product of an orthogonal matrix Q and an upper‑triangular matrix R. It’s one of the most versatile decompositions in linear algebra, especially for solving least‑squares problems and eigenvalue computations.

Breaks a matrix A (square or rectangular) into an orthogonal matrix and an upper‑triangular matrix:

A=QR

• Q → orthogonal (columns are orthonormal) (or unitary) matrix, meaning Q^TQ=I. Its columns form an orthonormal basis.
• R → upper‑triangular

If A is square and invertible, the decomposition is unique if we require the diagonal entries of R to be positive.

Several algorithms can produce the QR decomposition:

• Gram–Schmidt process → classical method, but numerically unstable.
• Householder reflections → stable and widely used in practice.
• Givens rotations → efficient for sparse matrices or when zeroing specific entries.

Used in least‑squares problems, numerical stability, and orthogonalization.

import numpy as np

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

Q, R = np.linalg.qr(A)

print("Q:\n", Q)
print("R:\n", R)

Applications

• Least Squares Problems: Instead of solving A^TAx=A^Tb, QR avoids squaring the condition number, making it numerically stable.
• Eigenvalue Computation: The QR algorithm iteratively applies QR decomposition to converge to eigenvalues.
• Orthonormal Basis Construction: Columns of Q provide an orthonormal basis for the column space of A.
• Numerical Stability: Preferred over LU decomposition in least‑squares contexts.

Cholesky Decomposition

Cholesky decomposition is a matrix factorization technique that expresses a symmetric, positive‑definite matrix as the product of a lower triangular matrix and its transpose. It is widely used because it is faster and more numerically stable than LU decomposition for such matrices.

For symmetric positive‑definite matrices:

• L → lower triangular matrix with strictly positive diagonal entries
• L^T → transpose of L

In the complex case (Hermitian matrices), the decomposition is:

where L^∗ is the conjugate transpose.

• Very fast and stable
• Used in optimization, Gaussian processes, covariance matrices

Properties

• Existence: Every symmetric positive‑definite matrix has a unique Cholesky decomposition.
• Efficiency: Roughly twice as efficient as LU decomposition for solving linear systems.
• Uniqueness: The decomposition is unique if we require positive diagonal entries in L.
• Semidefinite case: If A is only positive semidefinite, a decomposition still exists but may not be unique.

Algorithm (Step‑by‑Step)

For a matrix A, the entries of L are computed as:

This builds L row by row.

L = np.linalg.cholesky(A)

Example in Python

import numpy as np

A = np.array([[4, 12, -16],
              [12, 37, -43],
              [-16, -43, 98]])

L = np.linalg.cholesky(A)

print("L:\n", L)
print("Reconstructed A:\n", L @ L.T)

Applications

• Numerical linear algebra: Solving systems Ax=b efficiently.
• Optimization: Used in interior‑point methods and quadratic programming.
• Statistics: Sampling from multivariate normal distributions (covariance matrix factorization).
• Machine learning: Gaussian processes, Kalman filters, covariance estimation.
• Monte Carlo simulations: Efficiently generating correlated random variables.

Limitations

• Only works for symmetric positive‑definite matrices.
• If the matrix is not positive‑definite, the algorithm fails (square root of negative number).
• For semidefinite matrices, decomposition exists but may require pivoting and is not unique.

LU Decomposition (Not in NumPy, but in SciPy)

LU decomposition (or LU factorization) breaks a square matrix into the product of a lower triangular matrix L and an upper triangular matrix U. It is a fundamental tool in linear algebra for solving systems of equations, computing determinants, and finding matrix inverses.

For a square matrix A:

A=L⋅U

• L → lower triangular matrix (often with 1s on the diagonal).
• U → upper triangular matrix.

Sometimes a permutation matrix P is included to handle pivoting:

PA=LU

LU decomposition is essentially Gaussian elimination written in matrix form:

1. Start with matrix A.
2. Apply elimination steps to zero out entries below the diagonal.

• The multipliers used in elimination become entries of L.
• The resulting upper triangular matrix is U.

3. Pivoting (row swaps) may be required for numerical stability, leading to PA=LU.

Example

NumPy does not include LU decomposition directly. If you need it:

from scipy.linalg import lu
P, L, U = lu(A)

import numpy as np
from scipy.linalg import lu

A = np.array([[4, 3],
              [6, 3]])

P, L, U = lu(A)

print("P:\n", P)
print("L:\n", L)
print("U:\n", U)

Applications

Solving linear systems: Instead of solving Ax=b directly, factor A=LU.

• Solve Ly=b (forward substitution).
• Solve Ux=y (back substitution).

Matrix inverse: Compute inverse efficiently using LU decomposition.

Determinant: det⁡(A)=det⁡(L)⋅det⁡(U). Since det⁡(L)=1, determinant = product of diagonal entries of U.

Numerical methods: Used in iterative algorithms and simulations.

Variants

• Doolittle’s method: L has 1s on the diagonal.
• Crout’s method: U has 1s on the diagonal.
• Cholesky decomposition: Special case for symmetric positive‑definite matrices (A=LLT).

Limitations

• LU decomposition requires the matrix to be non‑singular (invertible).
• Pivoting may be necessary for stability (partial or complete pivoting).
• Not all matrices can be decomposed without pivoting.

Matrix Inverse and Pseudoinverse

These aren’t decompositions themselves, but they rely on them internally.

Inverse (square matrices)

A_inv = np.linalg.inv(A)

Pseudoinverse (any matrix, uses SVD)

A_pinv = np.linalg.pinv(A)

Solving Linear Systems

NumPy uses decomposition internally (LU or QR depending on the case):

x = np.linalg.solve(A, b)

polynomial mathematics

NumPy has a whole set of tools for polynomial mathematics under the numpy.polynomial module. These let you create, manipulate, and evaluate polynomials in a clean, object‑oriented way. Let’s go step by step.

1. Representing Polynomials

There are two main ways:

Using np.poly1d

import numpy as np

p = np.poly1d([1, -3, 2])   # x² - 3x + 2
print(p)
#    2
# 1 x - 3 x + 2

• Coefficients are given in descending powers.
• You can evaluate directly: p(5) → 12.

Using numpy.polynomial classes

NumPy provides specialized classes:

• Polynomial (coefficients in ascending order)
• Chebyshev, Legendre, Laguerre, Hermite (special orthogonal polynomials)

Example:

from numpy.polynomial import Polynomial

p = Polynomial([2, -3, 1])  # 2 - 3x + x²
print(p)
# 2.0 - 3.0 x + 1.0 x**2

build a polynomial directly from its roots

numpy.polynomial gives you a very clean way to build a polynomial directly from its roots using the function:

numpy.polynomial.Polynomial.fromroots()

This constructs the polynomial

with coefficients in ascending order (constant term first).

from numpy.polynomial import Polynomial

roots = [1, 2, 3]
p = Polynomial.fromroots(roots)

print(p)
# -6.0 + 11.0 x - 6.0 x**2 + 1.0 x**3

This matches the expansion of

Get the Coefficients

print(p.coef)
# [-6. 11. -6.  1.]

Using np.poly1d Instead

If you prefer the classic NumPy style:

import numpy as np

roots = [1, 2, 3]
p = np.poly1d(np.poly(roots))

print(p)
#    3     2
# 1 x - 6 x + 11 x - 6

polyval

polyval is NumPy’s classic function for evaluating a polynomial at specific points. It works with the traditional NumPy polynomial format, where coefficients are given in descending powers.

It evaluates a polynomial:

at one or many values of x.

import numpy as np

# Polynomial: 3x² + 2x + 1
coeffs = [3, 2, 1]

print(np.polyval(coeffs, 5))
# 86

Because:

Evaluating at Multiple Points

x = np.array([0, 1, 2, 3])
y = np.polyval([3, 2, 1], x)

print(y)
# [ 1  6 17 34]

Using with np.poly1d

p = np.poly1d([3, 2, 1])
print(np.polyval(p, 4))
# 57

• Coefficients must be in descending order.
• Works with scalars, lists, NumPy arrays.
• For more advanced polynomial math, prefer numpy.polynomial.Polynomial, which uses ascending‑order coefficients and has better numerical stability.

2. Polynomial Operations

You can add, subtract, multiply, and divide polynomials like normal objects.

p1 = np.poly1d([1, 2])   # x + 2
p2 = np.poly1d([1, -1])  # x - 1

print(p1 * p2)           # (x+2)(x-1) = x² + x - 2

3. Roots, Derivatives, Integrals

Roots:

p = np.poly1d([1, -3, 2])   # x² - 3x + 2
print(p.roots)              # [2, 1]

Derivative:

print(p.deriv())            # 2x - 3

Integral:

print(p.integ())            # (1/3)x³ - (3/2)x² + 2x

4. Fitting Data with Polynomials

NumPy can fit polynomials to data using least squares:

x = np.linspace(-1, 1, 10)
y = 2*x**2 + 3*x + 1 + np.random.randn(10)*0.1

coeffs = np.polyfit(x, y, deg=2)   # quadratic fit
p = np.poly1d(coeffs)

print("Fitted polynomial:", p)

np.polyfit is NumPy’s function for polynomial curve fitting. It finds the polynomial of a given degree that best fits your data in the least‑squares sense.

numpy.polyfit(x, y, deg)

• x → 1‑D array of x‑coordinates (independent variable)
• y → 1‑D array of y‑coordinates (dependent variable)
• deg → degree of the polynomial you want to fit

Returns:

• A 1‑D array of polynomial coefficients (highest power first).

import numpy as np
import matplotlib.pyplot as plt

# Sample data
x = np.linspace(-5, 5, 20)
y = 2*x**2 + 3*x + 1 + np.random.randn(20)*5   # noisy quadratic

# Fit a degree-2 polynomial
coeffs = np.polyfit(x, y, deg=2)
p = np.poly1d(coeffs)

print("Coefficients:", coeffs)
print("Fitted polynomial:\n", p)

# Plot
plt.scatter(x, y, label="Data")
plt.plot(x, p(x), color="red", label="Fitted curve")
plt.legend()
plt.show()

• np.polyfit uses least squares minimization.
• For high‑degree polynomials, results may be unstable (Runge’s phenomenon).
• If you want more numerically stable fitting, use numpy.polynomial.Polynomial.fit, which rescales the domain.

Applications

• Regression analysis: fitting curves to experimental data
• Trend detection: approximating nonlinear relationships
• Signal processing: smoothing noisy data
• Machine learning: feature engineering with polynomial terms

5. Special Polynomial Families

NumPy supports orthogonal polynomials useful in approximation and numerical analysis:

• numpy.polynomial.Chebyshev
• numpy.polynomial.Legendre
• numpy.polynomial.Laguerre
• numpy.polynomial.Hermite

Example:

from numpy.polynomial import Chebyshev

c = Chebyshev([0, 1])   # T1(x) = x
print(c(0.5))           # Evaluate at x=0.5