Mathematical Operations and Functions

In this chapter, we'll dive deep into NumPy's mathematical operations and functions, including basic arithmetic, trigonometric functions, statistical functions, linear algebra, and more.

Basic Arithmetic Operations

Element-wise Operations

import numpy as np

# Create sample arrays
array1 = np.array([1, 2, 3, 4, 5])
array2 = np.array([10, 20, 30, 40, 50])

print("Array 1:", array1)
print("Array 2:", array2)
print()

# Basic arithmetic operations
print("Addition array1 + array2:", array1 + array2)
print("Subtraction array2 - array1:", array2 - array1)
print("Multiplication array1 * array2:", array1 * array2)
print("Division array2 / array1:", array2 / array1)
print("Integer division array2 // array1:", array2 // array1)
print("Modulo array2 % array1:", array2 % array1)
print("Power array1 ** 2:", array1 ** 2)
print()

# Operations with scalar
scalar = 10
print(f"Array plus scalar array1 + {scalar}:", array1 + scalar)
print(f"Array times scalar array1 * {scalar}:", array1 * scalar)
print(f"Array divided by scalar array1 / {scalar}:", array1 / scalar)
print(f"Scalar minus array {scalar} - array1:", scalar - array1)

Mathematical Functions

import numpy as np

# Create test array
array = np.array([1, 4, 9, 16, 25])
print("Original array:", array)
print()

# Basic math functions
print("Square root np.sqrt():", np.sqrt(array))
print("Cube root np.cbrt():", np.cbrt(array))
print("Square np.square():", np.square(array))
print("Absolute value np.abs():", np.abs(array - 10))
print()

# Exponential and logarithmic functions
small_array = np.array([1, 2, 3, 4, 5])
print("Small array:", small_array)
print("Natural exponential np.exp():", np.exp(small_array))
print("Base-2 exponential np.exp2():", np.exp2(small_array))
print("Natural logarithm np.log():", np.log(small_array))
print("Base-10 logarithm np.log10():", np.log10(small_array))
print("Base-2 logarithm np.log2():", np.log2(small_array))
print()

# Rounding functions
float_array = np.array([1.2, 2.7, -1.5, -2.8, 3.0])
print("Float array:", float_array)
print("Floor np.floor():", np.floor(float_array))
print("Ceiling np.ceil():", np.ceil(float_array))
print("Round np.round():", np.round(float_array))
print("Truncate np.trunc():", np.trunc(float_array))

Trigonometric Functions

Basic Trigonometric Functions

import numpy as np

# Create angle array (radians)
angles_rad = np.array([0, np.pi/6, np.pi/4, np.pi/3, np.pi/2, np.pi])
print("Angles (radians):", angles_rad)
print("Angles (degrees):", np.degrees(angles_rad))
print()

# Basic trig functions
print("Sine np.sin():", np.sin(angles_rad))
print("Cosine np.cos():", np.cos(angles_rad))
print("Tangent np.tan():", np.tan(angles_rad))
print()

# Inverse trigonometric functions
values = np.array([0, 0.5, 0.707, 0.866, 1.0])
print("Values:", values)
print("Arcsine np.arcsin():", np.arcsin(values))
print("Arccosine np.arccos():", np.arccos(values))
print("Arctangent np.arctan():", np.arctan(values))
print()

# Hyperbolic functions
print("Hyperbolic sine np.sinh():", np.sinh([0, 1, 2]))
print("Hyperbolic cosine np.cosh():", np.cosh([0, 1, 2]))
print("Hyperbolic tangent np.tanh():", np.tanh([0, 1, 2]))

Angle Conversion and Special Functions

import numpy as np

# Angle and radian conversion
degrees = np.array([0, 30, 45, 60, 90, 180])
radians = np.radians(degrees)

print("Degrees:", degrees)
print("Converted to radians:", radians)
print("Converted back to degrees:", np.degrees(radians))
print()

# Special trigonometric functions
x = np.array([1, 2, 3])
y = np.array([1, 1, 4])

print("x:", x)
print("y:", y)
print("atan2(y, x):", np.arctan2(y, x))  # Quadrant-aware arctangent
print("hypot(x, y):", np.hypot(x, y))    # Calculate sqrt(x²+y²)
print()

# Angle normalization
angles = np.array([0, 90, 180, 270, 360, 450, -90])
print("Original angles:", angles)
print("Normalized to [0, 360):", angles % 360)
print("Normalized to [-180, 180):", ((angles + 180) % 360) - 180)

Statistical Functions

Basic Statistics

import numpy as np

# Create test data
np.random.seed(42)
data = np.random.normal(50, 10, 20)  # Normal distribution, mean 50, std 10
print("Test data:")
print(data)
print()

# Basic statistics
print(f"Mean np.mean(): {np.mean(data):.2f}")
print(f"Median np.median(): {np.median(data):.2f}")
print(f"Standard deviation np.std(): {np.std(data):.2f}")
print(f"Variance np.var(): {np.var(data):.2f}")
print(f"Minimum np.min(): {np.min(data):.2f}")
print(f"Maximum np.max(): {np.max(data):.2f}")
print(f"Range np.ptp(): {np.ptp(data):.2f}")  # Peak to peak
print(f"Sum np.sum(): {np.sum(data):.2f}")
print(f"Product np.prod(): {np.prod(data):.2e}")
print()

# Percentiles
percentiles = [25, 50, 75, 90, 95]
print("Percentiles:")
for p in percentiles:
    value = np.percentile(data, p)
    print(f"  {p}th percentile: {value:.2f}")

Multidimensional Array Statistics

import numpy as np

# Create 2D array
array_2d = np.random.randint(1, 10, (4, 5))
print("2D array:")
print(array_2d)
print()

# Global statistics
print(f"Global mean: {np.mean(array_2d):.2f}")
print(f"Global std: {np.std(array_2d):.2f}")
print()

# Statistics along axes
print("Row-wise statistics (axis=1):")
print(f"  Mean per row: {np.mean(array_2d, axis=1)}")
print(f"  Max per row: {np.max(array_2d, axis=1)}")
print(f"  Min per row: {np.min(array_2d, axis=1)}")
print()

print("Column-wise statistics (axis=0):")
print(f"  Mean per column: {np.mean(array_2d, axis=0)}")
print(f"  Max per column: {np.max(array_2d, axis=0)}")
print(f"  Min per column: {np.min(array_2d, axis=0)}")
print()

# Cumulative statistics
print("Cumulative sum np.cumsum():")
print(np.cumsum(array_2d, axis=1))  # Row-wise cumulative
print()

print("Cumulative product np.cumprod():")
print(np.cumprod(array_2d, axis=0))  # Column-wise cumulative

Correlation and Covariance

import numpy as np

# Create correlated data
np.random.seed(42)
x = np.random.normal(0, 1, 100)
y = 2 * x + np.random.normal(0, 0.5, 100)  # y correlated with x
z = np.random.normal(0, 1, 100)  # z uncorrelated with x, y

data = np.column_stack([x, y, z])
print(f"Data shape: {data.shape}")
print("First 5 rows:")
print(data[:5])
print()

# Correlation coefficient matrix
corr_matrix = np.corrcoef(data.T)
print("Correlation coefficient matrix:")
print(corr_matrix)
print()

# Covariance matrix
cov_matrix = np.cov(data.T)
print("Covariance matrix:")
print(cov_matrix)
print()

# Correlation between two variables
corr_xy = np.corrcoef(x, y)[0, 1]
print(f"Correlation between x and y: {corr_xy:.3f}")
print(f"Correlation between x and z: {np.corrcoef(x, z)[0, 1]:.3f}")
print(f"Correlation between y and z: {np.corrcoef(y, z)[0, 1]:.3f}")

Linear Algebra

Basic Matrix Operations

import numpy as np

# Create matrices
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
B = np.array([[9, 8, 7], [6, 5, 4], [3, 2, 1]])
v = np.array([1, 2, 3])

print("Matrix A:")
print(A)
print("Matrix B:")
print(B)
print("Vector v:", v)
print()

# Matrix multiplication
print("Matrix multiplication A @ B:")
print(A @ B)
print()

print("Matrix multiplication np.dot(A, B):")
print(np.dot(A, B))
print()

# Matrix-vector multiplication
print("Matrix-vector multiplication A @ v:")
print(A @ v)
print()

# Transpose
print("Transpose of A:")
print(A.T)
print()

# Trace (sum of diagonal elements)
print(f"Trace of A: {np.trace(A)}")
print(f"Trace of B: {np.trace(B)}")

Advanced Linear Algebra Operations

import numpy as np

# Create invertible matrix
A = np.array([[2, 1, 1], [1, 3, 2], [1, 0, 0]])
print("Matrix A:")
print(A)
print()

# Determinant
det_A = np.linalg.det(A)
print(f"Determinant det(A): {det_A:.3f}")
print()

# Matrix inverse
if det_A != 0:
    A_inv = np.linalg.inv(A)
    print("Inverse of A:")
    print(A_inv)
    print()
    
    # Verify A * A^(-1) = I
    identity = A @ A_inv
    print("A @ A^(-1) (should be identity matrix):")
    print(identity)
    print()

# Eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(A)
print("Eigenvalues:")
print(eigenvalues)
print("Eigenvectors:")
print(eigenvectors)
print()

# Matrix rank
rank_A = np.linalg.matrix_rank(A)
print(f"Matrix rank: {rank_A}")
print()

# Singular value decomposition
U, s, Vt = np.linalg.svd(A)
print("Singular value decomposition:")
print(f"U shape: {U.shape}")
print(f"Singular values: {s}")
print(f"Vt shape: {Vt.shape}")

Solving Linear Systems

import numpy as np

# Solve linear system Ax = b
# Example: 2x + y + z = 8
#          x + 3y + 2z = 20
#          x + 0y + 0z = 2

A = np.array([[2, 1, 1], [1, 3, 2], [1, 0, 0]])
b = np.array([8, 20, 2])

print("Coefficient matrix A:")
print(A)
print("Constant vector b:", b)
print()

# Solve using solve
x = np.linalg.solve(A, b)
print("Solution vector x:", x)
print()

# Verify solution
verification = A @ x
print("Verification A @ x:", verification)
print("Original b:", b)
print("Error:", np.abs(verification - b))
print()

# Least squares solution (for overdetermined systems)
A_over = np.array([[2, 1, 1], [1, 3, 2], [1, 0, 0], [0, 1, 1]])
b_over = np.array([8, 20, 2, 5])

print("Overdetermined system:")
print("A_over:")
print(A_over)
print("b_over:", b_over)

# Least squares solution
x_lstsq = np.linalg.lstsq(A_over, b_over, rcond=None)[0]
print("Least squares solution:", x_lstsq)

# Calculate residual
residual = A_over @ x_lstsq - b_over
print("Residual:", residual)
print(f"Sum of squared residuals: {np.sum(residual**2):.6f}")

Complex Number Operations

Complex Basics

import numpy as np

# Create complex array
complex_array = np.array([1+2j, 3-4j, 5+0j, 0+6j])
print("Complex array:", complex_array)
print("Data type:", complex_array.dtype)
print()

# Complex properties
print("Real part np.real():", np.real(complex_array))
print("Imaginary part np.imag():", np.imag(complex_array))
print("Magnitude np.abs():", np.abs(complex_array))
print("Phase angle np.angle():", np.angle(complex_array))
print("Conjugate np.conj():", np.conj(complex_array))
print()

# Complex operations
z1 = 3 + 4j
z2 = 1 - 2j

print(f"z1 = {z1}")
print(f"z2 = {z2}")
print(f"z1 + z2 = {z1 + z2}")
print(f"z1 * z2 = {z1 * z2}")
print(f"z1 / z2 = {z1 / z2}")
print(f"|z1| = {abs(z1)}")
print(f"z1* = {np.conj(z1)}")

Broadcasting

Understanding Broadcasting

import numpy as np

# Operations with different shaped arrays
array_2d = np.array([[1, 2, 3], [4, 5, 6]])
array_1d = np.array([10, 20, 30])
scalar = 100

print("2D array:")
print(array_2d)
print(f"Shape: {array_2d.shape}")
print()

print("1D array:", array_1d)
print(f"Shape: {array_1d.shape}")
print()

print(f"Scalar: {scalar}")
print()

# Broadcasting operations
print("2D array + 1D array:")
result_2d_1d = array_2d + array_1d
print(result_2d_1d)
print(f"Result shape: {result_2d_1d.shape}")
print()

print("2D array + scalar:")
result_2d_scalar = array_2d + scalar
print(result_2d_scalar)
print(f"Result shape: {result_2d_scalar.shape}")
print()

# Column vector broadcasting
column_vector = np.array([[1], [10]])
print("Column vector:")
print(column_vector)
print(f"Shape: {column_vector.shape}")

print("2D array * column vector:")
result_col = array_2d * column_vector
print(result_col)
print(f"Result shape: {result_col.shape}")

Practical Examples

Example 1: Signal Processing

import numpy as np

# Generate signal
t = np.linspace(0, 1, 1000)  # Time axis
frequency1 = 5  # 5Hz
frequency2 = 10  # 10Hz

# Combined signal
signal = np.sin(2 * np.pi * frequency1 * t) + 0.5 * np.sin(2 * np.pi * frequency2 * t)
noise = 0.2 * np.random.normal(0, 1, len(t))
noisy_signal = signal + noise

print(f"Signal length: {len(signal)}")
print(f"Signal statistics:")
print(f"  Mean: {np.mean(noisy_signal):.4f}")
print(f"  Standard deviation: {np.std(noisy_signal):.4f}")
print(f"  Maximum: {np.max(noisy_signal):.4f}")
print(f"  Minimum: {np.min(noisy_signal):.4f}")
print()

# Simple filtering (moving average)
window_size = 10
filtered_signal = np.convolve(noisy_signal, np.ones(window_size)/window_size, mode='same')

print("After filtering statistics:")
print(f"  Standard deviation: {np.std(filtered_signal):.4f}")
print(f"  SNR improvement: {np.std(noisy_signal)/np.std(filtered_signal):.2f}x")
print()

# Frequency domain analysis (simple power spectrum)
fft_result = np.fft.fft(signal)
power_spectrum = np.abs(fft_result)**2
frequencies = np.fft.fftfreq(len(signal), t[1] - t[0])

# Find main frequency components
positive_freq_idx = frequencies > 0
positive_freqs = frequencies[positive_freq_idx]
positive_power = power_spectrum[positive_freq_idx]

# Find frequency with maximum power
max_power_idx = np.argmax(positive_power)
dominant_frequency = positive_freqs[max_power_idx]

print(f"Dominant frequency component: {dominant_frequency:.1f} Hz")

Example 2: Financial Data Analysis

import numpy as np

# Simulate stock price data
np.random.seed(42)
days = 252  # Trading days per year
initial_price = 100
daily_returns = np.random.normal(0.001, 0.02, days)  # Daily returns
prices = initial_price * np.cumprod(1 + daily_returns)

print(f"Stock price analysis ({days} trading days):")
print(f"  Initial price: ${initial_price:.2f}")
print(f"  Final price: ${prices[-1]:.2f}")
print(f"  Total return: {(prices[-1]/initial_price - 1)*100:.2f}%")
print()

# Return statistics
print("Daily return statistics:")
print(f"  Average daily return: {np.mean(daily_returns)*100:.3f}%")
print(f"  Return std: {np.std(daily_returns)*100:.3f}%")
print(f"  Largest gain: {np.max(daily_returns)*100:.2f}%")
print(f"  Largest loss: {np.min(daily_returns)*100:.2f}%")
print()

# Risk metrics
# 1. Annualized volatility
annual_volatility = np.std(daily_returns) * np.sqrt(252)
print(f"Annualized volatility: {annual_volatility*100:.2f}%")

# 2. Maximum drawdown
running_max = np.maximum.accumulate(prices)
drawdown = (prices - running_max) / running_max
max_drawdown = np.min(drawdown)
print(f"Maximum drawdown: {max_drawdown*100:.2f}%")

# 3. Sharpe ratio (assuming 2% risk-free rate)
risk_free_rate = 0.02
excess_returns = daily_returns - risk_free_rate/252
sharpe_ratio = np.mean(excess_returns) / np.std(excess_returns) * np.sqrt(252)
print(f"Sharpe ratio: {sharpe_ratio:.3f}")
print()

# Moving averages
window_short = 20  # 20-day MA
window_long = 50   # 50-day MA

ma_short = np.convolve(prices, np.ones(window_short)/window_short, mode='valid')
ma_long = np.convolve(prices, np.ones(window_long)/window_long, mode='valid')

print(f"Technical indicators:")
print(f"  Current price: ${prices[-1]:.2f}")
print(f"  20-day MA: ${ma_short[-1]:.2f}")
print(f"  50-day MA: ${ma_long[-1]:.2f}")

# Simple trading signal
if len(ma_short) > 0 and len(ma_long) > 0:
    if ma_short[-1] > ma_long[-1]:
        signal = "Buy signal"
    else:
        signal = "Sell signal"
    print(f"  Trading signal: {signal}")

Performance Optimization

Vectorization vs Loops

import numpy as np
import time

# Compare vectorized and loop performance
size = 1000000
array1 = np.random.random(size)
array2 = np.random.random(size)

# Using loop
start_time = time.time()
result_loop = np.zeros(size)
for i in range(size):
    result_loop[i] = array1[i] * array2[i] + np.sin(array1[i])
loop_time = time.time() - start_time

# Using vectorization
start_time = time.time()
result_vectorized = array1 * array2 + np.sin(array1)
vectorized_time = time.time() - start_time

print(f"Array size: {size:,}")
print(f"Loop method time: {loop_time:.4f} seconds")
print(f"Vectorized method time: {vectorized_time:.4f} seconds")
print(f"Performance improvement: {loop_time/vectorized_time:.1f}x")
print(f"Results equal: {np.allclose(result_loop, result_vectorized)}")

Chapter Summary

In this chapter, we learned:

• Basic arithmetic operations and mathematical functions
• Trigonometric and inverse trigonometric functions
• Statistical functions and descriptive statistics
• Linear algebra operations
• Complex number operations
• Broadcasting mechanism
• Practical examples (signal processing, financial analysis)
• Performance optimization techniques

Next Steps

In the next chapter, we'll learn NumPy's advanced features, including random number generation, polynomials, Fourier transforms, and more.

Exercises

1. Create a function to calculate standardized scores (z-score) of an array
2. Implement a simple moving average filter
3. Write a function to calculate cosine similarity between two vectors
4. Use linear algebra to solve a multiple linear regression problem
5. Implement a simple image convolution function