111+ Mathematics Problems with Python

for AI and Quantitative Finance

1 Computational Levels
The problems in this book are organized into three computational levels, allowing you to
progressively develop your programming and numerical computing skills while solving the
same mathematical problems. The goal is to solve all the problems from Level 1 to Level 3
progressively. The different computational levels are as follows.

Level 1: Python Loop Thinking

Level 1 introduces the foundations of computational problem solving using core Python. You
learn to express algorithms through explicit control flow using loops, conditionals, functions,
and basic data structures such as lists, dictionaries, and tuples. At this stage, computa-
tion is performed one element at a time, emphasizing algorithmic reasoning and step-by-step
execution.

Level 2: Python Collection Thinking

Level 2 develops a higher level of abstraction by treating collections of data as objects to
be transformed rather than iterated over manually. You learn to write cleaner and more
expressive programs using list comprehensions, nested comprehensions, zip, and enumerate.
The emphasis shifts from explicit iteration to applying transformations over entire collections
while remaining within the Python language.

Level 3: Tensor Vectorized Thinking

Level 3 introduces vectorized numerical computing using NumPy. Instead of writing explicit
loops, you learn to express computations directly on arrays using vectorized operations, broad-
casting, masking, reshaping, axis-based reductions, and matrix operations. The focus is on
writing mathematical expressions that operate on entire arrays simultaneously, providing the
foundation for scientific computing, machine learning, artificial intelligence, & quant finance.
The goal is to perform batched computation as if a list of inputs are there instead of one.
2 Crash Course: Three Levels of Computation
Computation begins with Level 1: Python Loop Thinking, where all operations are
performed explicitly one element at a time. At this stage, data is treated as simple scalars
stored inside lists, and computation is expressed using for loops, while loops, and basic
control flow such as if/else. The mental model is strictly sequential: each value is processed
individually, and results are accumulated manually using variables, dictionaries, or lists. This
level builds fundamental algorithmic thinking, but it is slow and does not scale well for
numerical computing.
The next stage is Level 2: Python Collection Thinking, where the focus shifts from
individual elements to entire collections. Instead of manually iterating step-by-step, we begin
transforming datasets using list comprehensions, zip, enumerate, and functional patterns
like map and filter. Computation is now expressed as operations over structures rather
than loops over elements. The mental model becomes: “think in transformations, not
iterations.” This level significantly improves code clarity and prepares the transition to
mathematical array reasoning, but computation is still fundamentally Python-based and not
yet optimized for numerical performance.
Finally, Level 3: Tensor Vectorized Thinking introduces true numerical computing
through ndarray objects. Here, computation is expressed directly on arrays without explicit
loops, using broadcasting, elementwise operations, axis-based reductions, reshaping, masking,
and matrix algebra. Functions such as mean, sum, and argmax operate over entire dimensions
simultaneously, allowing computations to be executed in optimized compiled code. The men-
tal model becomes: “write mathematics directly on arrays.” This is the foundation
of scientific computing, machine learning, and quantitative finance, where performance and
abstraction depend on vectorization rather than iteration.

Level 1 (Loops) → Level 2 (Collections) → Level 3 (Vectorized Arrays)

The fundamental shift across these levels is not syntax, but abstraction: computation evolves
from element-wise reasoning, to structure-wise transformation, and finally to vectorized math-
ematical execution. This progression is the core intellectual bridge between Python program-
ming and modern scientific computing.
 44+ Linear Algebra Problems with Python
for AI and Quantitative Finance

1 Easy Problems
1.1 Matrix-Vector Multiplication.
Matrix-vector multiplication forms the core computational kernel of linear layers in deep
learning neural networks and asset allocation weight mappings. Given a transformation matrix
A ∈ Rm×n and a column vector x ∈ Rn , we seek to calculate the transformed vector y ∈
Rm where each element yi is computed as Pnthe dot product of the i-th row of A with x.
Mathematically, this is expressed as yi = j=1 Aij xj . Input: A nested list of real numbers
representing an m × n matrix A, and a list of real numbers representing a vector x of length
n. Output: A list of real numbers representing the resulting vector y of length m.

Examples of Different Computational Levels

def mat_ ve c_m ul _l evel1 (A , x ) : # LEVEL 1

m = len ( A )
n = len ( x )
y = [0.0] * m
for i in range ( m ) :
total = 0.0
for j in range ( n ) :
total += A [ i ][ j ] * x [ j ]
y [ i ] = total
return y

def mat_ ve c_m ul _l evel2 (A , x ) : # LEVEL 2

return [ sum ( a_ij * x_j for a_ij , x_j in zip ( row , x ) ) for row in
A]
 import numpy as np

def mat_ ve c_m ul _l evel3 (A , x ) : # LEVEL 3 ( Batched Computation )

A = np . array ( A ) # Shape : (B , M , N ) or (M , N )
x = np . array ( x ) # Shape : (B , N , 1) or (N ,)

# Performing batched matrix multiplication over batch dimension

y = np . matmul (A , x )
return y

# Example input with a batch of 2 matrices and 2 vectors

A_batch = [[[1 , 2] , [3 , 4]] , [[2 , 0] , [1 , 3]]]
x_batch = [[1 , 2] , [3 , 4]]
y = mat_ ve c_ mu l_ le vel3 ( A_batch , x_batch )

1.2 Vector Dot Product.

The dot product represents the foundational projection of one directional coordinate vector
onto another, acting as the fundamental metric for calculating unscaled cosine similarities or
expected portfolio returns. Given two vectors u and v of identical length n, the dot product
maps their P
elements pairwise to a single scalar landscape. Mathematically, it is defined as
n
d = u · v = i=1 ui vi . Input: Two lists of real numbers u = [u1 , . . . , un ] and v = [v1 , . . . , vn ].
Output: A real scalar number representing the inner product d ∈ R.

1.3 Vector Norms (L1 and L2 ).

Vector norms provide a rigorous definition of distance and size within multi-dimensional vec-
tor spaces, forming the structural basis for Lasso (L1 ) and Ridge (L2 ) regularization penalties
in machine learning models. Given a spatial vector x ∈ Rn , the L1 norm evaluates the ab-
solute Manhattan distance, Pwhile the L2 norm evaluates
pPn the straight-line Euclidean distance.
n 2
Mathematically, ∥x∥1 = i=1 |x i | and ∥x∥ 2 = i=1 xi . Input: A list of real numbers
X = [x1 , . . . , xn ]. Output: A tuple of two real numbers representing the L1 norm and L2
norm.

1.4 Matrix Transposition.

Transposition flips a matrix over its diagonal, switching its row and column indices. This util-
ity is critical for restructuring dimension alignments during operations like backpropagation
or updating historical financial correlation data. Given a matrix A ∈ Rm×n , the transpose
matrix AT ∈ Rn×m satisfies the element-wise coordinate mapping constraint ATji = Aij . In-
put: A nested list of real numbers representing an m × n matrix A. Output: A nested list of
real numbers representing the transformed n × m matrix AT .
1.5 Matrix Trace Calculation.
The trace extracts geometric structural properties from square transformation arrays, serving
as an invariant scalar property that equals the sum of the matrix’s eigenvalues. Given a
square matrix A ∈ Rn×n , the trace Tr(A) sums the entries along its primary Pnmain diagonal
from top-left to bottom-right. Mathematically, this is formalised as Tr(A) = i=1 Aii . Input:
A nested list of real numbers representing a square n × n matrix A. Output: A real scalar
number representing the computed trace Tr(A).

1.6 Cosine Similarity Metric.

Cosine similarity measures the angular alignment between two directional vectors independent
of their scale, making it a key metric for text embeddings and item filtering algorithms.
Given two non-zero vectors u, v ∈ Rn , the similarity score evaluates the cosine of the angle
between them, bounding the result within [−1, 1]. Mathematically, sim(u, v) = ∥u∥u·v 2 ∥v∥2
=
P
√P 2u√
i vi
. Input: Two lists of real numbers uandv of equal length. Output: A real number
vi2
P
ui
representing the calculated cosine similarity score.

1.7 Frobenius Norm of a Matrix.

The Frobenius norm quantifies the total size or energy of an entry matrix by mapping its
multi-dimensional elements down to a single scalar distance, which helps measure matrix-
level model weight regularization errors. Given an underlying matrix A ∈ Rm×n , it calculates
the square root
qPof thePsum of all squared coordinate entries. Mathematically, this is expressed
m n 2
as ∥A∥F = i=1 j=1 Aij . Input: A nested list of real numbers representing an m × n
matrix A. Output: A real scalar number representing the calculated Frobenius norm.

1.8 Vector Outer Product.

The outer product combines two spatial vectors to construct a rank-one matrix structural
map, which forms the building block for constructing covariance updates and cross-attention
alignment steps in transformer layers. Given a vector u ∈ Rm and a vector v ∈ Rn , the outer
product matrix M = u ⊗ v ∈ Rm×n contains elements defined by the entry-wise product
Mij = ui vj . Input: A list of real numbers u of length m, and a list of real numbers v of length
n. Output: A nested list of real numbers representing the m × n outer product matrix.

1.9 Matrix Symmetry Verification.

Symmetric structures are common in physics and optimization, ensuring that covariance and
asset correlation matrices have real eigenvalues and orthogonal eigenvectors. Given a square
matrix A ∈ Rn×n , we evaluate whether it matches its own transpose across a small numerical
tolerance threshold ϵ. Mathematically, the matrix is symmetric if it satisfies the constraint
|Aij − Aji | < ϵ for all index pairs i, j. Input: A nested list representing a square matrix A
and a real threshold ϵ > 0. Output: A boolean value (True or False) indicating whether the
matrix is symmetric.
1.10 Hadamard Element-wise Product.
The Hadamard product performs element-wise multiplication of matching array cells, which is
a key operation for modulating features through gating mechanisms in networks like LSTMs
or GRUs. Given two matrices A, B ∈ Rm×n with identical dimensions, the resulting matrix
C = A ⊙ B ∈ Rm×n is constructed by multiplying elements with matching indices: Cij =
Aij · Bij . Input: Two nested lists of real numbers representing m × n matrices A and B.
Output: A nested list of real numbers representing the element-wise product matrix C.

1.11 Linear Combination of Vectors.

Linear combinations form the foundation of vector spaces, allowing us to build composite
spaces and calculate weighted portfolio asset positions by scaling and combining individual
base vectors. Given a set array sequence of vectors v1 , v2 , . . . , vk ∈ Rn along with a corre-
Pk set of scaling scalars c1 , c2 , . . . , ck ∈ R, we calculate the structural target vector
sponding
w = i=1 ci vi . Input: A nested list of real numbers representing k vectors of length n, and a
list of k scalar values. Output: A list of real numbers representing the combined target vector
w.

1.12 Matrix Diagonal Component Extraction.

Extracting or manipulating the primary diagonal of a matrix is a common preprocessing step
for initializing iterative solvers or isolating variance entries from empirical covariance matrices.
Given a square input matrix A ∈ Rn×n , the goal is to isolate its main diagonal elements into
a single 1D vector d ∈ Rn where di = Aii for i ∈ {1, . . . , n}. Input: A nested list of real
numbers representing an n × n square matrix A. Output: A list of real numbers representing
the isolated diagonal entries.

1.13 Scalar Multiplication of a Matrix.

Scaling structural properties uniformly maps to uniform adjustments of portfolio exposures
or tuning the learning rate factors inside deep weight systems. Given a matrix A ∈ Rm×n and
a real scalar factor c ∈ R, we scale each element independently. Mathematically, the scaled
output matrix components satisfy Bij = c · Aij . Input: A nested list representing an m × n
matrix A and a real number scalar c. Output: A nested list representing the scaled matrix B.

1.14 Vector Addition and Subtraction.

Combining independent movement channels underpins translation steps and baseline updates
across linear optimization models. Given two coordinate arrays u, v ∈ Rn of matching length,
we combine components across single index fields. Mathematically, the resulting spatial array
maps out element-wise sum wi = ui ± vi . Input: Two lists of real numbers u and v of identical
length. Output: A tuple containing two lists representing the structural sum and difference
vectors.
1.15 Standard Basis Vector Generator.
Generating elementary directional units allows coordinate systems to map discrete feature
coordinates into absolute vector spaces. Given a targeting vector space dimension size n
and an active component anchor index k, we construct a unit coordinate vector ek ∈ Rn .
Mathematically, ek = [x1 , . . . , xn ]T where xi = 1 if i = k else 0. Input: An integer tracking
total dimensional depth n ≥ 1 and an integer targeting index k (0 ≤ k < n). Output: A list
of real numbers representing the standard basis vector ek .

2 Medium Problems
2.1 Determinant of a 3 × 3 Matrix.
The determinant measures the volume scaling factor of a linear transformation, indicating
whether a matrix mapping is invertible or if it collapses the vector space into a lower di-
mension. Given a square matrix A ∈ R3×3 , we calculate its scalar determinant via cofac-
tor expansion along its primary row boundary. Mathematically, det(A) = A11 (A22 A33 −
A23 A32 ) − A12 (A21 A33 − A23 A31 ) + A13 (A21 A32 − A22 A31 ). Input: A nested list of real num-
bers representing a 3 × 3 matrix A. Output: A real scalar number representing the computed
determinant det(A).

2.2 Matrix-Matrix Multiplication.

Matrix multiplication represents the composition of sequential linear transformations, serving
as the foundational computational bottleneck for training deep neural architectures and pro-
cessing multi-asset risk projections. Given two conformal matrices A ∈ Rm×n and B ∈ Rn×p ,
m×p
the resulting matrix C = AB ∈ RP is computed by taking the dot product of the rows of
n
A with the columns of B: Cik = j=1 Aij Bjk . Input: A nested list representing an m × n
matrix A, and a second nested list representing an n × p matrix B. Output: A nested list
representing the resulting m × p matrix C.

2.3 Solving Linear Systems via Gaussian Elimination.

Solving systems of linear equations in the form Ax = b is a fundamental problem in linear
algebra, used to balance financial portfolios, solve ordinary differential equations, and compute
network traffic flows. Given an invertible coefficient matrix A ∈ Rn×n and a target vector
b ∈ Rn , this algorithm uses forward elimination to reduce the augmented matrix [A|b] to upper
triangular form, followed by back-substitution to isolate the coordinate solutions. Input: A
nested list representing matrix A, and a list representing vector b. Output: A list of real
numbers representing the solution vector x that satisfies Ax = b.

2.4 Gram-Schmidt Orthogonalization Process.

The Gram-Schmidt process transforms a set of linearly independent vectors into an orthogonal
or orthonormal basis that spans the exact same geometric subspace, which is a key step
for stabilizing numerical algorithms and performing QR factorizations. Given a sequence
of vectors v1 , . . . , vk ∈ RnP
, we iteratively subtract projections onto previously computed
i−1 v ·u
reference axes: ui = vi − j=1 ∥ui j ∥j2 uj , then normalize each vector to unit length. Input:
A nested list of real numbers representing a collection of k independent vectors of length n.
Output: A nested list of real numbers representing the resulting orthonormalized vectors.

2.5 Eigenvalues of a Symmetric 2 × 2 Matrix.

Eigenvalues identify the scaling factors along a transformation’s principal axes of distor-
tion, making them essential for principal component analysis
 and verifying the stability of
a b
dynamical systems. Given a symmetric matrix A = , we find its eigenvalues by solv-
b c
ing the characteristic equation det(A − λI) = 0, which simplifies to the quadratic equation
λ2 − Tr(A)λ + det(A) = 0. Input: A nested list of real numbers representing a symmetric 2 × 2
matrix A. Output: A tuple of two real numbers representing the sorted eigenvalues λ1 ≤ λ2 .

2.6 Singular Value Decomposition (SVD) for Dimensionality Reduc-

tion.
Singular Value Decomposition factorizes an arbitrary data matrix into principal singular struc-
tural components, serving as the mathematical foundation for latent semantic analysis, image
compression, and low-rank matrix approximations. Given a matrix A ∈ Rm×n , SVD fac-
torizes it into A = U ΣV T , where U and V are orthogonal matrices containing the left and
right singular vectors, and Σ is a diagonal matrix containing the singular values sorted in
descending order. Input: A nested list representing an m × n matrix A, and an integer k
specifying the target low-rank truncation limit. Output: A tuple of three arrays representing
the truncated components Uk , Σk , and VkT .

2.7 Principal Component Analysis (PCA) Projection.

PCA projects high-dimensional datasets onto their directions of maximum variance, filter-
ing out noise and reducing dimensionality for downstream classification and pricing models.
Given a data matrix X ∈ Rn×d , the process centers the data by subtracting the mean vec-
1 T
tor, computes the empirical covariance matrix Σ = n−1 Xcentered Xcentered , extracts its top
eigenvectors, and projects the centered data onto this lower-dimensional subspace. Input: A
nested list representing an n×d data matrix X, and an integer k specifying the target number
of dimensions (k < d). Output: A nested list representing the projected n × k data matrix.

2.8 Matrix Quadratic Form and Definiteness Verification.

Quadratic forms map a vector to a scalar value through a central structural weight matrix,
which is a key framework for evaluating utility functions, portfolio risk variances (xT Σx), and
second-order optimization P gradients. Given a symmetric matrix A ∈ Rn×n , we evaluate the
T
scalar value q = x Ax = i,j Aij xi xj and analyze its definiteness by checking the signs of
its eigenvalues or principal minors. Input: A nested list representing a square matrix A, and
a list representing a coordinate vector x. "Positive Definite", "Indefinite").
2.9 Vector Projection onto a Subspace Axis.
Vector projection breaks down a data vector into orthogonal components relative to a reference
directional axis, forming the foundational computational kernel for least-squares approxima-
tions and hyperplanar classification boundaries. Given a target vector v ∈ Rn and a reference
direction vector u ∈ Rn , the parallel projection vector p minimizes the distance ∥v − p∥2 and
v·u
is calculated using the formula p = ∥u∥ 2 u. Input: A list representing the target vector v,
2
and a second list representing the reference basis vector u. Output: A list of real numbers
representing the projected vector coordinate p.

2.10 Linear Independence and Matrix Rank Estimation.

Estimating matrix rank identifies the number of linearly independent rows or columns in a
matrix, which determines whether a linear system has a unique solution or if a data matrix
contains redundant, collinear features. Given an arbitrary matrix A ∈ Rm×n , we transform it
to row echelon form using Gaussian elimination or compute its singular values to determine
the rank, which equals the number of non-zero rows or singular values above a small numerical
tolerance ϵ. Input: A nested list representing an m × n matrix A, and a small real tolerance
threshold ϵ > 0. Output: An integer representing the computed rank of the matrix.

2.11 Least Squares Solution via Normal Equations.

When a linear system Ax = b has more equations than variables, it is overdetermined and
typically has no exact solution. The method of least squares finds an approximate solution
x̂ that minimizes the squared residual error ∥Ax − b∥22 , which is the standard optimization
criterion for linear regression models. This optimal solution satisfies the normal equations,
given by AT Ax̂ = AT b, which can be solved directly as x̂ = (AT A)−1 AT b. Input: A nested
list representing an m × n design matrix A, and a list representing an m-dimensional target
vector b. Output: A list of real numbers representing the least-squares parameter solution
vector x̂.

2.12 Moore-Penrose Pseudoinverse Calculation.

The Moore-Penrose pseudoinverse generalizes the concept of a matrix inverse to non-square or
rank-deficient matrices, providing a unified framework for solving general linear systems and
finding least-squares solutions. Given a matrix A ∈ Rm×n with singular value decomposition
A = U ΣV T , its pseudoinverse A+ ∈ Rn×m is computed as A+ = V Σ+ U T , where Σ+ is
constructed by taking the reciprocal of each non-zero diagonal singular value and transposing
the matrix. Input: A nested list of real numbers representing an arbitrary m × n matrix A.
Output: A nested list of real numbers representing the computed pseudoinverse matrix A+ .

2.13 Mahalanobis Distance Metric.

The Mahalanobis distance measures the distance between a data point and a distribution in
a multi-dimensional space, accounting for correlations between variables to serve as a robust
metric for anomaly detection and pattern recognition. Given a data vector x ∈ Rd , a mean
distribution focus vector µ ∈ Rd , and the dataset’s covariance matrix
p Σ∈R
d×d
, the distance
metric is calculated using the quadratic form equation DM (x) = (x − µ) Σ−1 (x − µ). In-
2

put: A list representing vector x, a list representing vector µ, and a nested list representing
the positive-definite covariance matrix Σ. Output: A real scalar number representing the
calculated Mahalanobis distance.

2.14 Matrix Exponential (eA ) via Power Series.

Matrix exponentials are essential for solving systems of linear differential equations and mod-
eling continuous-time Markov transition processes in financial asset pricing networks. Given
a square matrix A ∈ Rn×n , the matrix exponential eA is defined by a convergent power series
expansion
P∞ analogous to the scalar exponential function. Mathematically, this is expressed as
1 k 1 2 1 3
eA = k=0 k! A = I + A + 2! A + 3! A + . . . , which can be approximated numerically by
truncating the series at a sufficiently high order K. Input: A nested list representing a square
matrix A, and an integer K specifying the number of expansion terms to compute. Output:
A nested list of real numbers representing the approximated exponential matrix eA .

2.15 Householder Reflection Transformation.

Householder reflections use mirror-image transformations across a hyperplanar boundary to
zero out specific elements of a vector, serving as the foundational computational block for
stable QR factorizations and tridiagonalizing matrices. Given an initial vector x ∈ Rn , we
x−∥x∥2 e1
construct a unit householder vector v = ∥x−∥x∥ 2 e1 ∥2
to build an orthogonal reflection matrix
H = I − 2vv T , which maps x onto a multiple of the standard basis vector e1 . Input: A
list of real numbers representing a spatial coordinate vector x. Output: A nested list of real
numbers representing the square orthogonal Householder matrix H.

2.16 LU Decomposition with Partial Pivoting.

LU decomposition factors a square matrix into the product of a lower triangular matrix L
and an upper triangular matrix U . This decomposition simplifies downstream computations,
allowing linear systems to be solved efficiently via forward and backward substitution. To
ensure numerical stability and prevent division by zero, rows are swapped during elimination
using a permutation matrix P , yielding the standard decomposed form P A = LU . Input: A
nested list of real numbers representing a square n × n matrix A. Output: A tuple of three
nested lists representing the permutation matrix P , the unit lower triangular matrix L, and
the upper triangular matrix U .

2.17 Cholesky Decomposition for Sampling.

Decomposing positive-definite matrices isolates directional variance structures, enabling cor-
related Monte Carlo paths or generating multivariate Gaussian samples. Given a symmetric
positive-definite covariance matrix Σ ∈ Rn×n , we seek a lower triangular matrix
q L such that
Pi−1
Σ = LLT . Mathematically, the diagonal entries are computed as Lii = Σii − k=1 L2ik .
Input: A nested list representing a square symmetric positive-definite matrix Σ. Output: A
nested list representing the unique lower triangular matrix factor L.

2.18 Matrix Condition Number Calculation.

The condition number gauges the sensitivity of linear equations to measurement errors, in-
dicating how numerical instabilities propagate during matrix inversion or backpropagation
steps. Given an invertible array A, the metric tracks the product of structural norms between
the base and inverse shapes. Mathematically, κ(A) = ∥A∥ · ∥A−1 ∥, often evaluated using the
ratio of extreme singular values σσmax
min
. Input: A nested list representing a square matrix A.
Output: A real scalar number representing the matrix condition number κ(A).

2.19 QR Decomposition via Givens Rotations.

QR decomposition factorizes a matrix into an orthogonal operator Q and an upper trian-
gular array R, which stabilizes least-squares evaluations. Givens rotations achieve this fac-
torization by systematically zeroing out sub-diagonal entries using localized coordinate ro-
tations. Mathematically, a plane rotation matrix G(i, j, θ) modifies rows i and j such that
Gii = c, Gjj = c, Gij = s, Gji = −s, where c = cos θ and s = sin θ. Input: A nested list
of real numbers representing a square n × n matrix A. Output: A tuple of two nested lists
representing the orthogonal matrix Q and the upper triangular matrix R.

2.20 Power Iteration for Dominant Eigenpair.

Power iteration isolates the single largest eigenvalue and its corresponding eigenvector, which
provides the mathematical foundation for ranking algorithms like PageRank and tracking
dominant risk drivers. Given a square matrix A ∈ Rn×n and a random starting vector
x0 , the algorithm iteratively multiplies the vector by the matrix and normalizes it: xk+1 =
Axk xT
k Axk
∥Axk ∥2 . The Rayleigh quotient xT k xk
then converges to the dominant eigenvalue λmax . Input:
A nested list representing a square matrix A and an integer tracking maximum iterations
N . Output: A tuple containing the real dominant eigenvalue λ and a list representing its
eigenvector x.

2.21 Gerschgorin Disc Theorem Eigenvalue Bounds.

The Gerschgorin Disc Theorem provides quick analytical boundaries for bounding the regional
landscape of eigenvalues without running costly iterative diagonalization solvers. Given a
complex square matrix A ∈ Rn×n , every eigenvalue λ lies within at least one closed disc
D(Aii , Ri ) in the complex plane. Mathematically, the center of each disc is the diagonal
entry Aii , and
P its radius is the sum of the absolute values of the remaining entries in that
row: Ri = j̸=i |Aij |. Input: A nested list of real numbers representing a square matrix
A. Output: A list of tuples, where each tuple contains the real center value Aii and the
corresponding real radius Ri .
2.22 Kronecker Product of Two Matrices.
The Kronecker product interleaves two arbitrary matrices to construct a larger block ma-
trix, which maps out complex interaction systems and structures multi-asset tensor states
in quantitative finance models. Given matrices A ∈ Rm×n and B ∈ Rp×q , the Kronecker
product C = A ⊗ B ∈ Rmp×nq is formed by scaling the entire matrix B by each individual
element of A. Mathematically, this is expressed as a block matrix where block (i, j) is given
by Aij B. Input: Two nested lists representing matrix A and matrix B. Output: A nested
list representing the large tensor product matrix C.

2.23 Matrix Null Space Basis Extraction.

The null space isolates the subspace of non-trivial input vectors that map directly to zero,
which helps detect linear dependencies and identify redundant dimensions in portfolio alloca-
tion frameworks. Given a matrix A ∈ Rm×n , we seek the solution space containing all vectors
x that satisfy the constraint Ax = 0. This basis can be extracted by identifying the free
variables after reducing the matrix to reduced row echelon form (RREF). Input: A nested list
representing an m×n matrix A. Output: A nested list of vectors representing an orthonormal
basis for the kernel null space.

2.24 Vector Cross Product (3D Space).

The cross product calculates a vector that is perpendicular to two input vectors in three-
dimensional space, providing a foundational metric for computing geometric torque, structural
normals, and angular momentum paths. Given vectors u, v ∈ R3 , the cross product vector
w = u × v ∈ R3 is computed using the determinant of a formal matrix. Mathematically,
w = [u2 v3 − u3 v2 , u3 v1 − u1 v3 , u1 v2 − u2 v1 ]T . Input: Two lists of real numbers u and v of
length 3. Output: A list of real numbers representing the perpendicular vector w.

2.25 Vectorization of a Matrix (vec operator).

The vectorization operator flattens a multi-dimensional matrix into a single column vector
by stacking its columns sequentially, which converts matrix equations into equivalent vector
forms for optimization solvers. Given an arbitrary matrix A ∈ Rm×n , the vectorized column
vector v = vec(A) ∈ Rmn maps each coordinate entry using the index transform formula
vi+(j−1)m = Aij . Input: A nested list of real numbers representing an m × n matrix A.
Output: A single flat list of real numbers representing the column-stacked vector v.

2.26 Projection Matrix Construction for a Subspace.

A projection matrix maps any vector in a larger space onto its closest point within a specific
subspace, which forms the computational basis for orthogonal regressions and dimension
reduction filters. Given a matrix X ∈ Rm×n whose columns span the target subspace, the
orthogonal projection matrix P is symmetric and idempotent (P 2 = P ). Mathematically, it
is constructed using the formula P = X(X T X)−1 X T . Input: A nested list representing a
matrix X whose columns span the target subspace. Output: A square nested list representing
the orthogonal projection matrix P .
2.27 Sherman-Morrison Formula Update.
The Sherman-Morrison formula computes the inverse of a matrix after a rank-one update
without running a costly full inversion from scratch, which speeds up online learning al-
gorithms and recursive portfolio optimization. Given an invertible matrix A ∈ Rn×n and
column vectors u, v ∈ Rn , the inverse of the updated matrix (A + uv T ) is computed as
−1 T −1
(A + uv T )−1 = A−1 − A1+vuv A
T A−1 u . Input: A nested list representing the original inverse
−1
matrix A , and two lists representing vectors u and v. Output: A nested list representing
the updated inverse matrix.

2.28 Symmetric Rank-1 Update (SR1).

Symmetric rank-one updates adjust approximation matrices sequentially while preserving
symmetry, which forms a key step in quasi-Newton optimization techniques like the SR1
algorithm for updating Hessian approximations. Given a symmetric matrix Bk , a change in
gradient yk , and a step vector sk , the updated matrix Bk+1 must satisfy the secant equation
T
Bk+1 sk = yk . Mathematically, it is updated as Bk+1 = Bk + (yk −B k sk )(yk −Bk sk )
(yk −Bk sk )T sk
. Input: A
symmetric matrix Bk , a step vector sk , and a displacement vector yk . Output: A nested list
representing the updated symmetric matrix Bk+1 .

2.29 Vandermonde Matrix Construction.

Vandermonde matrices map a vector of points to a polynomial coordinate space, which is
a key step for performing polynomial interpolation, least-squares curve fitting, and signal
processing transformations. Given a vector of geometric points x = [x1 , . . . , xn ], the resulting
matrix V ∈ Rn×m raises each input point to sequential powers. Mathematically, the matrix
entries are defined by the power relation Vij = xj−1
i . Input: A list of real numbers x of length
n, and an integer specifying the target number of column degrees m. Output: A nested list
representing the structured n × m Vandermonde matrix.

2.30 Rayleigh Quotient Bounds Estimation.

The Rayleigh quotient measures how much a matrix scales a vector across different directions,
and its minimum and maximum values bound the matrix’s entire eigenvalue spectrum. Given
a symmetric matrix A ∈ Rn×n and a non-zero vector x ∈ Rn , the Rayleigh quotient R(A, x)
maps the vector to a scalar value. Mathematically, it is calculated using the ratio R(A, x) =
xT Ax
xT x
. Input: A square symmetric matrix A and a non-zero coordinate evaluation vector x.
Output: A real scalar number representing the calculated Rayleigh quotient value R(A, x).

2.31 Matrix Commutator Calculation.

The commutator measures the degree of non-commutativity between two linear transforma-
tions, which determines whether two transformations can be applied in any order or if their
sequence alters the final state. Given two square matrices A, B ∈ Rn×n , the commutator op-
erator [A, B] evaluates the difference between their forward and reverse compositions. Math-
ematically, it is defined as [A, B] = AB − BA. Input: Two nested lists representing square
matrices A and B of equal dimensions. Output: A nested list representing the resulting matrix
commutator.
 44+ Probability Problems with Python
for AI and Quantitative Finance

1 Easy Problems
1.1 Mean, Median, Mode.
Descriptive statistics form the core foundation of data analysis, summarizing central tenden-
cies of datasets in AI feature engineering and finance. Given an empirical list of P numbers
n
representing observed samples, we seek to calculate the arithmetic mean µ = n1 i=1 xi ,
the median which splits the sorted sample array into equal halves, and the mode defined
as the most frequently occurring value arg maxk freq(xk ). Input: A list of real numbers
X = [x1 , x2 , . . . , xn ]. Output: A tuple of three real numbers representing the mean µ ∈ R,
median m ∈ R, and mode M ∈ R.

Examples of Different Computational Levels

def me a n _ m e d i a n _ m o de _l e ve l 1 ( X ) : # LEVEL 1
n = len ( X )
total = 0
for x in X :
total += x
mean = total / n # mean
sorted_X = sorted ( X )
if n % 2 == 1:
median = sorted_X [ n // 2]
else :
median = ( sorted_X [ n // 2 - 1] + sorted_X [ n // 2]) / 2 #
median
freq = {}
for x in X :
if x in freq :
freq [ x ] += 1
else :
freq [ x ] = 1
mode = max ( freq , key = freq . get ) # mode
return mean , median , mode
 from collections import Counter
def me a n _ m e d i a n _ m o de _l e ve l 2 ( X ) : # LEVEL 2
n = len ( X )
mean = sum ( X ) / n # Mean
X_sorted = sorted ( X ) # Median
mid = n // 2
median = ( X_sorted [ mid ] if n % 2 else ( X_sorted [ mid - 1] +
X_sorted [ mid ]) / 2)
mode = Counter ( X ) . most_common (1) [0][0] # Mode
return mean , median , mode

import numpy as np

def me a n _ m e d i a n _ m o de _l e ve l 3 ( X ) : # LEVEL 3 ( Batched Computation )

X = np . array ( X )
mean = np . mean (X , axis =1) # mean
median = np . median (X , axis =1) # median
modes = [] # mode ( no vectorized operation exist so for loop is
used )
for row in X :
values , counts = np . unique ( row , return_counts = True )
modes . append ( values [ np . argmax ( counts ) ])
mode = np . array ( modes )
return mean , median , mode
X = [[1 ,2 ,2 ,3 ,4] ,
[5 ,6 ,7 ,8 ,9] ,
[2 ,2 ,3 ,3 ,4]]
mean , median , mode = me an _ me d ia n_ m od e _l ev e l3 ( X )

1.2 Sample Variance & Standard Deviation.

Quantifying volatility or data dispersion is essential in risk management and standardizing
machine learning model features. Given a list of numerical data points, the sample variance s2
measures the average squared deviation from the sample mean x̄ using Bessel’s correction n−1,
while the standard deviation s returns this P measure to the original unit
√ scale. Mathematically,
1 n
these parameters are defined as s2 = n−1 (x
i=1 i − x̄) 2
and s = s2 . Input: A list of real
numbers X = [x1 , x2 , . . . , xn ]. Output: A tuple of two real numbers representing the sample
variance s2 ∈ R and sample standard deviation s ∈ R.

1.3 Percentiles / Quantiles.

In quantitative finance, Value at Risk (VaR) models depend explicitly on mapping data
cutoffs using percentiles to isolate tail risk. Given an empirical array of values and a target
percentile q, the objective is to sort the data and apply a linear interpolation scheme to
calculate the score v below which q% of the values fall. Mathematically, after sorting the
array such that x(1) ≤ x(2) ≤ · · · ≤ x(n) , the percentile rank map maps a continuous value
q
v = x(i) + (r − i)(x(i+1) − x(i) ) where r = 1 + 100 (n − 1) and i = ⌊r⌋. Input: A list of real
numbers X = [x1 , . . . , xn ] and a real number percentile q ∈ [0, 100]. Output: A real number
representing the calculated percentile value v ∈ R.

1.4 Expected Value (Discrete Distribution).

The expected value acts as the foundational long-run average outcome for random variables,
mapping directly to asset pricing in financial markets. Given a finite set of outcomes along
with their verified individual discrete probabilities, the expectation is formalised as the dot
product ofPnthe outcome vector and probability vector. Mathematically,
Pn this is expressed as
E[X] = i=1 x i · P (X = x i ), subject to the constraint i=1 P (X = xi ) = 1. Input: A
list of unique real outcomes X = [x1 , . . . , xn ] and a corresponding list of real probabilities
P = [p1 , . . . , pn ]. Output: A real number representing the expected value E[X] ∈ R.

1.5 Bernoulli PMF & Moments.

The Bernoulli distribution forms the binary atomic building block for modeling random fail-
ure/success gates, classification nodes, or defaults in credit risk. Given a success probability
parameter p, we evaluate a binary outcome x ∈ {0, 1} and determine the structural moments.
Mathematically, the PMF is f (x; p) = px (1 − p)1−x , while the foundational moments are
defined cleanly by the expected value E[X] = p and the variance Var(X) = p(1 − p). Input:
A real number success probability p ∈ [0, 1] and an integer target state x ∈ {0, 1}. Output: A
tuple containing a real PMF probability, a real mean, and a real variance.

1.6 Binomial PMF.

The binomial distribution models discrete processes counting independent, repeated trials,
such as predicting total option exercise frequencies under a static tree. Given the count
of trials n, success rate p, and target success number k, we evaluate the static probability
at a single discrete point using
combinatorics. Mathematically, this is computed using the
PMF formula f (k; n, p) = nk pk (1 − p)n−k , where the binomial coefficient is calculated as
n n!


k = k!(n−k)! . Input: An integer number of trials n, an integer number of successes k, and
a real probability p ∈ [0, 1]. Output: A real number representing the exact point probability
P (X = k) ∈ [0, 1].

1.7 Geometric PMF & Mean.

The geometric distribution tracks the trial index of the first isolated success in a sequence of
independent Bernoulli trials, mimicking search procedures or default timing. Given a single
success parameter p, we compute the point probability of finding the first success on trial
k along with its theoretical asymptotic average. Mathematically, the PMF is derived as
f (k; p) = (1 − p)k−1 p and the expected mean is given by E[X] = p1 . Input: A real number
success probability p ∈ (0, 1] and an integer trial index k ≥ 1. Output: A tuple containing a
real PMF value and a real expected mean value.
1.8 Expected Value and Variance of an n-Sided Die.
Evaluating symmetric multi-faced discrete uniform systems forms the primary basis for build-
ing unbiased randomized gaming engines or simple non-parametric baseline simulations. Given
a balanced n-sided die containing uniform faces from 1 to n, we compute the long-run expec-
tation and spread analytically. Mathematically, the expected value is E[X] = n+12 , and the
n2 −1
exact analytical variance simplifies directly to Var(X) = 12 . Input: An integer represent-
ing total faces n ≥ 1. Output: A tuple of two real numbers representing the expected value
E[X] ∈ R and variance Var(X) ∈ R.

1.9 Sampling Distribution of the Mean.

The structural framework of statistical inference relies heavily on predicting how sample
aggregates behave relative to their underlying parent parameter populations. Given a known
population mean µ, population standard deviation σ, and a set draw size n, we evaluate the
parameters of the sampling distribution of the sample mean X̄. Mathematically, according
to foundational sampling theory, the expected mean remains E[X̄] = µ, while the standard
error is σx̄ = √σn . Input: A real mean µ, a real standard deviation σ, and an integer sample
size n ≥ 1. Output: A tuple of two real numbers representing the expected mean and the
standard error.

1.10 Calculate Conditional Probability from Data.

Empirical conditional probabilities allow AI algorithms to construct native frequency baseline
matrices using raw tabular categorical observations. Given a historical data log of parallel
pairs containing binary status reports for events A and B, we count sub-frequencies to find the
proportion of active outcomes. Mathematically,
Pn this evaluates the basic frequentist conditional
P (A∩B) I(Ai =1∧Bi =1)
division identity P (A|B) = P (B) = i=1
P n , where I is the indicator function.
i=1 I(Bi =1)
Input: A list of binary coordinate pairs [(A1 , B1 ), (A2 , B2 ), . . . , (An , Bn )] where Ai , Bi ∈
{0, 1}. Output: A real number representing the conditional empirical probability P (A|B) ∈
[0, 1].

1.11 Compute Posterior Probability using Bayes’ Theorem.

Bayes’ Theorem provides the primary mechanism for sequentially updating diagnostic belief
states when an agent encounters new evidence. Given a base prior probability P (A) along
with true positive likelihoods P (B|A) and false positive likelihoods P (B|Ac ), we compute the
posterior belief P (A|B). Mathematically, applying standard conditional probability yields
P (B|A)P (A)
P (A|B) = P (B|A)P (A)+P (B|Ac )(1−P (A)) . Input: A real prior probability P (A), a real likeli-
hood P (B|A), and a real likelihood P (B|Ac ), all bounded in [0, 1]. Output: A real number
representing the updated posterior probability P (A|B) ∈ [0, 1].

1.12 Simulate Two-Dice Sum Distribution.

Simulating discrete combinatoric outcome spaces provides a robust framework for validat-
ing exact theoretical distributions against empirical frequencies. Given a target simulation
count N , we randomly roll two distinct six-sided dice, calculate their sum, and tabulate the
normalized frequency distribution across all possible sums from 2 to P 12. Mathematically, let
N
D1,i , D2,i ∼ U{1, 6}, the objective is to approximate P (S = s) ≈ N1 i=1 I(D1,i + D2,i = s)
for s ∈ {2, 3, . . . , 12}. Input: An integer representing the total number of simulation iterations
N ≥ 1. Output: A dictionary mapping integer sums s ∈ {2, . . . , 12} to their real empirical
frequencies.

2 Medium Problems
2.1 Binomial Distribution Probability.
Evaluating localized or bounded success counts within complex operational pipelines helps
check quality bounds or systematic trade risk profiles. Given total discrete trials n, a single-
trial success likelihood p, and a specific exact target threshold k, we compute the explicit
point density function of this binomial structure. Mathematically, this is expressed as P (X =
k) = nk pk (1 − p)n−k , which requires managing large factorial limits securely or utilizing

stable log-gamma scaling transforms when parameters scale high. Input: An integer total
trial count n, an integer success target k, and a real probability p ∈ [0, 1]. Output: A real
number representing the calculated exact point probability.

2.2 Normal Distribution PDF Calculator.

The Gaussian continuous distribution is a cornerstone of modern financial engineering models
and deep learning optimization assumptions. Given a specific evaluation site x, along with
the underlying mean µ and standard deviation variance scale σ, we compute the structural
density coordinate. Mathematically, the normal
 probability
 density function is defined by the
1 (x−µ)2
analytical equation f (x; µ, σ) = σ√2π exp − 2σ2 . Input: A real number evaluation coor-
dinate x, a real number distribution mean µ, and a positive real number standard deviation
σ > 0. Output: A real number representing the continuous density value f (x) ∈ [0, ∞).

2.3 Poisson Distribution Probability Calculator.

The Poisson distribution is widely used to model the frequency of rare operational failures,
network packet bursts, or sudden credit defaults within a fixed time window. Given a contin-
uous arrival scale intensity rate λ and a target event count k, we compute the exact discrete
point density. Mathematically, the probability mass function is calculated using the formula
k −λ
P (X = k) = λ k! e
, requiring careful handling of exponential bounds to avoid numerical
overflow. Input: A real number intensity rate parameter λ > 0 and an integer event count
threshold k ≥ 0. Output: A real number representing the exact point probability.

2.4 Simulate Markov Chain Transitions.

Markov state transitions serve as basic frameworks for structural asset regime shifts, consumer
health migrations, or predictive text sequences. Given an explicit square stochastic transition
matrix P , a vector describing the initial state probabilities v, and a target time horizon
length N , we track individual state transitions over time. Mathematically, at step t, the next
state is sampled from the categorical
P distribution defined by row Xt of matrix P , satisfying
Pij = P (Xt+1 = j|Xt = i) and j Pij = 1. Input: A square matrix P ∈ RM ×M , an initial
probability distribution vector v ∈ RM , and an integer total steps parameter N ≥ 1. Output:
A list of N integers tracking the sequence of simulated state indices.

2.5 Calculate KL Divergence Between Two Multivariate Gaussian

Distributions.
Kullback-Leibler (KL) divergence provides information-geometric metrics for measuring the
divergence between a candidate model and a target distribution in continuous variational
autoencoders. Given the continuous vector means µ1 , µ2 and structural covariance matrices
Σ1 , Σ2 of two separate multivariate Normal systems, we calculate h this directional discrepancy.
|Σ2 |
1
Mathematically, the divergence is expressed as DKL (p1 ||p2 ) = 2 log |Σ 1|
− d + Tr(Σ−1
2 Σ1 ) + (µ2 −
where d is the vector dimension. Input: Two real vectors µ1 , µ2 ∈ Rd and two positive-definite
matrices Σ1 , Σ2 ∈ Rd×d . Output: A real number representing the calculated information di-
vergence DKL ∈ [0, ∞).

2.6 Chi-square Probability Distribution.

Chi-square structures are vital for executing continuous parameter sample testing and ver-
ifying goodness-of-fit across empirical models. Given the degrees of freedom parameter k
and an integration limit point x, we evaluate the continuous cumulative distribution function
(CDF). Mathematically, this is defined by the regularized lower incomplete gamma function,
R x/2 k/2−1 −t
expressed as F (x; k) = γ(k/2,x/2)
Γ(k/2)
1
= Γ(k/2) 0
t e dt, evaluated for non-negative do-
mains where x ≥ 0. Input: An integer degrees of freedom parameter k ≥ 1 and a non-negative
real calculation coordinate x ∈ [0, ∞). Output: A real number representing the integrated
cumulative probability F (x) ∈ [0, 1].

2.7 Conditional Probability from Joint Distribution.

Extracting individual conditional profiles from a discrete joint distribution allows us to analyze
dependencies between categorical variables. Given a two-dimensional matrix containing the
joint probability mass function (PMF) values for variables X and Y , along with specific
query indices x and y, we extract the slice profile. Mathematically, we evaluate the relation
P (X = x|Y = y) = P (X=x,Y =y)
P (Y =y) , where the marginal probability is computed by summing
across the row: P (Y = y) = x′ P (X = x′ , Y = y). Input: A joint probability mass matrix
P
M ∈ RR×C and target coordinate integers x and y. Output: A real number representing the
conditional probability P (X = x|Y = y) ∈ [0, 1].

2.8 Central Limit Theorem Simulation.

The Central Limit Theorem (CLT) shows why the normal distribution appears frequently
in nature: the sum of independent, identically distributed variables converges to a Gaus-
sian distribution. Given a sample draw size n and an execution iteration limit M , we re-
peatedly generate samples from a continuous uniform distribution, average them, and out-
put the empirical collection. Mathematically, for each iterationPj ∈ {1, . . . , M }, we sample
n
uj,1 , . . . , uj,n ∼ U(0, 1) and compute the sample mean x̄j = n1 i=1 uj,i . Input: An integer
tracking total generated iterations M and an integer tracking internal sample size size n.
Output: A list of M real numbers containing the empirical distribution of sample means.

2.9 Compute Covariance from Joint PMF.

Covariance measures the linear joint directional behavior of two random variables, which
is a key parameter for portfolio diversification and asset allocation. Given a discrete joint
PMF matrix mapping outcomes across two coordinate dimensions, we compute the expected
product of deviations
P P from their respective means. Mathematically, this is expressed as
Cov(X, Y ) = i j xP i yj P (X
P = xi , Y = yj )−E[X]E[Y ], where P
the marginal
P expectations are
computed as E[X] = i xi j P (X = xi , Y = yj ) and E[Y ] = j yj i P (X = xi , Y = yj ).
Input: A list of outcomes X, a list of outcomes Y , and a joint probability mass matrix
P ∈ R|X|×|Y | . Output: A real number representing the covariance Cov(X, Y ) ∈ R.

2.10 Compute Total Probability using Law of Total Probability.

The Law of Total Probability allows us to compute the global probability of an event by
partitioning the sample space into distinct, conditional sub-scenarios. Given a vector of
conditional probabilities P (A|Bi ) along with a vector of mutually exclusive prior probabilities
P (Bi ), we calculate the unconditional probability P (A).PMathematically, this is computed
m
using the total probability partition Psum formula P (A) = i=1 P (A|Bi )P (Bi ), subject to the
m
structural partitioning constraint i=1 P (Bi ) = 1. Input: A list of conditional probabilities
PA|B = [pA|B1 , . . . , pA|Bm ] and a list of structural priors PB = [pB1 , . . . , pBm ]. Output: A real
number representing the total integrated probability P (A) ∈ [0, 1].

2.11 Hypergeometric Distribution PMF.

The Hypergeometric distribution models random sampling without replacement, where each
selection alters the success probability of subsequent draws. Given a total population size N , a
total number of success states within that population K, a sample size n, and a target success
count k, we compute the exact probability. Mathematically, the probability mass function is
(K )(N −K )
given by P (X = k) = k Nn−k , evaluated within the constraints max(0, n − (N − K)) ≤ k ≤
(n)
min(n, K). Input: Four integers representing population size N , sub-population success count
K, sample count n, and target draw success count k. Output: A real number representing
the exact hypergeometric probability.

2.12 Birthday Problem Probability.

The birthday paradox illustrates how counterintuitive probability can be, with applications
ranging from hash collision limits in cryptography to security risk profiling. Given a group
size n inside a standard year, we compute the probability that at least two individuals share
the same calendar birthday. Mathematically, we compute the complement of the probability
 Qn−1
that all birthdays are unique: P (at least one match) = 1 − i=0 365−i
365 , for group sizes where
n ≤ 365. Input: An integer representing the total number of individuals in the group n ≥ 1.
Output: A real number representing the probability of at least one shared birthday.

2.13 Negative Binomial Distribution Probability.

The negative binomial distribution models the number of independent Bernoulli trials needed
to achieve a target number of successes, which is useful for predicting operational duration
or marketing thresholds. Given a target success count r, a total trial index milestone k,
and a success probability parameter p, we compute the probability that the r-th success
occurs exactly on the k-th

r trial. k−r
Mathematically, the probability mass function is calculated
as P (X = k) = k−1 r−1 p (1 − p) , valid for conditions where k ≥ r. Input: An integer
target success count r, an integer target milestone trial index k, and a real success probability
p ∈ [0, 1]. Output: A real number representing the point probability.

2.14 Monte Carlo Estimate of Expected Value in a Bidding Game.

When analytical solutions are difficult to derive, Monte Carlo simulations provide an effec-
tive way to estimate expectations in strategic bidding models or option contracts. Given a
simulation sample budget N and a static bid level, we simulate random opponent strate-
gies to computePN the average payoff. Mathematically, the expected value is approximated by
E[g(X)] ≈ N1 i=1 g(bid, Yi ), where Yi represents random opponent profiles drawn from a
specified distribution, and g is the payoff function. Input: An integer tracking simulation
loops N ≥ 1 and a real number player bid value b ∈ R. Output: A real number representing
the estimated empirical expected payoff.

2.15 Poisson PMF & CDF.

A comprehensive analysis of Poisson processes requires evaluating both individual point event
probabilities and cumulative risk thresholds over time. Given a continuous intensity arrival
scale rate λ and a target event count threshold k, we compute both the point density and
k −λ
the cumulative probability. Mathematically, the point probability is f (k; λ) = λ k!
e
, and the
Pk λi e−λ
cumulative distribution function is F (k; λ) = i=0 i! . Input: A real number intensity
parameter λ > 0 and an integer event count index k ≥ 0. Output: A tuple containing a real
PMF probability and a real CDF probability.

2.16 Chi-Square Test.

The Chi-Square independence test is a non-parametric method used to determine whether two
categorical variables are independent by comparing observed frequencies with expected base-
lines. Given a contingency matrix of observed counts, we compute the chi-square test statistic
PR PC (O −E )
to evaluate discrepancies. Mathematically, the statistic is calculated as χ2 = i=1 j=1 ijEij ij
P P
Oik m Omj
where the expected frequency matrix elements are defined as Eij = kP
Okm . Input: A
k,m

contingency matrix of observed counts O ∈ RR×C filled with non-negative numbers. Output:
A real number representing the computed test statistic χ2 ∈ [0, ∞).
2.17 One-Sample t-Test.
The one-sample t-test determines whether a sample mean significantly differs from a hypoth-
esized population mean, which is useful for evaluating trading strategy performance against
a benchmark. Given an empirical list of numbers and a baseline mean µ0 , we calculate the
t-statistic using the sample standard error. Mathematically, the test statistic is computed
as t = x̄−µ√ 0 , where x̄ is the sample mean, s is the sample standard deviation, and n is the
s/ n
sample size. Input: A list of empirical real numbers X = [x1 , . . . , xn ] and a real benchmark
mean µ0 ∈ R. Output: A real number representing the calculated test statistic t ∈ R.

2.18 Bootstrap Mean & Confidence Interval.

Bootstrapping is a non-parametric resampling technique used to estimate the variability of
a statistic and compute confidence intervals without making strong parametric assumptions.
Given a dataset and a resampling budget B, we draw samples with replacement to estimate
the lower and upper bounds of the mean’s
Pn confidence interval. Mathematically, we generate B
empirical bootstrap means x̄∗b = n1 i=1 x∗b,i , sort them, and extract the values at the α2 and
1 − α2 quantiles. Input: A list of numbers X = [x1 , . . . , xn ] and an integer tracking bootstrap
draw cycles B ≥ 1. Output: A tuple of two real numbers representing the lower and upper
confidence interval bounds.

2.19 Empirical Conditional Probability.

Empirical conditional probability measures the observed frequency of a target event given
that a conditioning event occurs, which helps assess dependencies in messy datasets. Given
an unformatted list of categorical pairs, we calculate the conditional probability by filtering
for thePconditioning event. Mathematically, this evaluates the empirical ratio P (A = 1|B =
n
I(A =1∧Bi =1)
1) = i=1 Pn i , returning zero if the denominator conditioning set is empty. Input:
i=1 I(Bi =1)
A list of paired attribute states [(A1 , B1 ), . . . , (An , Bn )] where Ai , Bi ∈ {0, 1}. Output: A real
number representing the empirical conditional probability P (A = 1|B = 1) ∈ [0, 1].

2.20 Discrete Random Variable Transformation.

Transforming random variables allows us to derive the distribution of a new variable, such
as mapping asset returns to squared volatility metrics in finance. Given a discrete variable
X defined by its probability mass function, we compute the probability distribution of a new
variable Y = g(X) = X 2 . Mathematically, the new probability massPfunction is calculated by
accumulating the probabilities of the original values: P (Y = y) = x:x2 =y P (X = x). Input:
A dictionary mapping real values xi to their respective probabilities P (X = xi ). Output: A
dictionary mapping transformed squared values yj to their calculated probabilities P (Y = yj ).

2.21 Empirical Variance of Conditional Distribution.

Conditional variance quantifies the residual uncertainty of a variable after accounting for a
conditioning factor, a key concept in autoregressive conditional heteroskedasticity (ARCH)
models. Given a dataset of paired observations (Xi , Yi ), we calculate the sample variance of
X withinP a specific slice where Y = y. Mathematically, this is expressed as Var(X|Y = y) =
1 2
|Iy |−1 i∈Iy (xi − x̄Iy ) , where Iy = {i : yi = y} and x̄Iy is the mean of that conditional subset.
Input: A list of paired data coordinates [(x1 , y1 ), . . . , (xn , yn )] and a target conditioning value
y. Output: A real number representing the empirical conditional variance Var(X|Y = y) ∈
[0, ∞).

2.22 Joint PMF Normalization Check.

Validating the structure of data matrices is an essential preprocessing step to ensure they
represent valid probability spaces before running downstream inference models. Given a two-
dimensional matrix of non-negative real numbers, we verify if the elements sum to 1 within a
small numerical tolerance. Mathematically, the matrix represents a valid joint distribution if
PR PC
Mij ≥ 0 for all i, j and satisfies the normalization constraint 1 − i=1 j=1 Mij < ϵ, where
ϵ is a small threshold. Input: A matrix of real numbers M ∈ RR×C and a small real tolerance
threshold ϵ > 0. Output: A boolean value (True or False) indicating whether the matrix is
a valid probability distribution.

2.23 Markov Chain n-step Transition Probability.

Long-term multi-period state predictions in Markov systems require projecting the transi-
tion matrix over multiple time horizons. Given an initial single-step transition matrix P
and a target step horizon n, we compute the n-step transition matrix, which represents the
probabilities of transitioning between states over n steps. Mathematically, according to the
Chapman-Kolmogorov equations, this matrix power sequence is computed via repeated matrix
multiplication: P (n) = P n = P × P × · · · × P . Input: A square transition matrix P ∈ RM ×M
| {z }
n times
and an integer multi-step horizon parameter n ≥ 1. Output: A square matrix representing
the n-step transition probabilities P n ∈ RM ×M .

2.24 Order Statistic: Minimum of Samples.

Order statistics are used to analyze extreme events, such as estimating the time-to-failure of
a system or modeling extreme risk in financial portfolios. Given an independent sampling
process, we run a simulation to estimate the expected value of the minimum value across a
sample of size n. Mathematically, given M simulation runs where each run draws n inde-
pendent values Xj,1 , . . . , Xj,n from a distribution F , we estimate the mean of the first order
1
PM
statistic: E[X(1) ] ≈ M j=1 min(Xj,1 , . . . , Xj,n ). Input: An integer tracking total simulation
cycles M ≥ 1 and an integer tracking the sample draw size n ≥ 1. Output: A real number
representing the estimated expected minimum value.

2.25 Order Statistic: Median Simulation.

Robust statistical models often use the sample median rather than the mean to minimize the
impact of outliers in noisy environments. Given a sample size n and a simulation budget M ,
we generate repeated samples from a base distribution to estimate the expected value of the
median. Mathematically, for each simulation run j ∈ {1, . . . , M }, we sort the n independent
draws to find the sample median mj = median(Xj,1 , . . . , Xj,n ), and approximate the expected
1
PM
value as E[Xmed ] ≈ M j=1 mj . Input: An integer tracking total simulation cycles M ≥ 1
and an integer tracking sample size n ≥ 1. Output: A real number representing the estimated
expected median value.

2.26 Law of Total Variance.

The law of total variance decomposes the total variance of a system into intra-group and
inter-group components, providing a framework for analyzing hierarchical structures. Given a
dataset containing paired observations of a target variable and its grouping labels, we compute
and compare the conditional variance components. Mathematically, we verify the variance
decomposition identity: Var(X) = E[Var(X|Y )] + Var(E[X|Y ]), processing the conditional
terms across the grouped subsets. Input: A list of paired observations [(x1 , y1 ), . . . , (xn , yn )]
where yi represents categorical grouping labels. Output: A tuple of two real numbers repre-
senting the total variance and the sum of the decomposed variance components.

2.27 Empirical PDF Estimation (Histogram).

Density estimation constructs continuous probability profiles from discrete empirical data
points without assuming an underlying parametric distribution. Given a dataset and a spec-
ified number of bins, we partition the data range and compute the normalized density for
each bin. Mathematically, for a bin spanning the interval [L, R) with width W = R − L, the
c
normalized empirical density is calculated as fbin = n·W , where c is the count of data points
falling within that bin and n is the total number of samples. Input: A list of empirical real
numbers X = [x1 , . . . , xn ] and an integer specifying the number of bins B ≥ 1. Output: A
list of B real numbers representing the normalized probability densities.

2.28 Continuous Random Variable Expectation Approximation.

When analytical integration is difficult, numerical approximations allow us to estimate the
expected values of continuous random variables. Given an intensity parameter λ and a sample
budget N , we use Monte Carlo integration to estimate the expected value of an exponential
distribution. Mathematically, we generate N independent uniform samples ui ∼ U(0, 1),
convert them using the inversePtransform sampling method xi = − λ1 log(1 − ui ), and compute
n
the sample mean: E[X] ≈ N1 i=1 xi . Input: An integer tracking total sample draws N ≥ 1
and a real rate parameter λ > 0. Output: A real number representing the approximated
continuous expectation.

2.29 Correlation Matrix Computation.

Correlation matrices summarize linear dependencies across multivariate datasets, serving as
a core input for portfolio optimization and dimensionality reduction models. Given an input
data matrix, we compute the Pearson correlation coefficient for all pairs of feature columns.
Cov(X ,X )
Mathematically, the elements of the correlation matrix R are calculated as Rij = σX σiX j ,
i j
where Cov represents the sample covariance and σ represents the sample standard deviation
of columns i and j. Input: A data matrix X ∈ Rn×d containing n samples and d features.
Output: A symmetric matrix representing the correlation coefficients R ∈ Rd×d , with diagonal
elements equal to 1.

3 Hard Problems
3.1 First Passage Probability Simulation.
First passage times help quantify structural risk by determining the probability that a stochas-
tic process crosses a critical barrier within a given time horizon. Given a Markov chain
transition matrix, a designated start state, and a target horizon T , we use simulation to
estimate the probability of hitting a specific state for the first time. Mathematically, we
simulate paths to estimate the probability P (τ ≤ T |X0 = sstart ), where the random variable
τ = inf{t ≥ 1 : Xt = starget } represents the first passage time index. Input: A transition
matrix P ∈ RM ×M , an integer start state, an integer target state, and an integer time horizon
limit T . Output: A real number representing the estimated first passage probability.

3.2 Bayesian Update with Multiple Evidence.

In complex environments like automated trading or robotics, algorithms must sequentially
update their internal beliefs as they receive a stream of noisy signals. Given an initial prior
probability distribution across a set of hypotheses and a sequence of observed signals, we
iteratively compute the updated posterior distribution. Mathematically, for a sequence of
independent pieces of evidence E = [e1 , . . . , ek ], the sequential update rule is formulated as
P (Hj ) k
Q
m=1 P (em |Hj )
P (Hj |e1 , . . . , ek ) = P Qk . Input: A vector of prior probabilities P (H) ∈
r P (Hr ) m=1 P (em |Hr )
M M ×V
R and a likelihood matrix L ∈ R mapping hypotheses to signal probabilities. Output:
A vector representing the final normalized posterior probability distribution P (H|E) ∈ RM .

3.3 Law of Large Numbers Convergence Speed.

The Law of Large Numbers guarantees long-term convergence, but analyzing the rate of
convergence is essential for understanding error bounds in numerical simulations. Given an
expanding sample size sequence up to a maximum limit N , we evaluate how quickly the
empirical sample mean converges to the true theoretical mean. Mathematically,Pn we
track
the structural convergence error as a function of the sample size: ϵ(n) = n1 i=1 xi − µ ,
and estimate the asymptotic error rate parameter α where ϵ(n) ≈ c · n−α . Input: An integer
tracking the maximum sample simulation limit N ≥ 100. Output: A real number representing
the calculated empirical convergence rate exponent α.

3.4 Joint Probability of Dependent Continuous Random Variables

via Monte Carlo Integration.
Evaluating multivariate probability bounds analytically becomes mathematically intractable
when random variables are bound by complex joint dependencies rather than simple indepen-
dence. Given two continuous random variables X and Y whose joint distribution is defined
by a dependent transformation of independent standard Gaussian variables, we use a simula-
tion framework to compute the joint probability that their values simultaneously fall within
a specific circular bound. Mathematically, given independent p variables Z1 , Z2 ∼ N (0, 1), we
construct the dependent system X = Z1 and Y = ρZ1 + 1 − ρ2 Z2 based on a correlation
parameter ρ, and approximate the probability P (X 2 + Y 2 ≤ R2 ) by evaluating the ratio of
PN
simulated coordinate pairs (xi , yi ) that satisfy the condition N1 i=1 I(x2i + yi2 ≤ R2 ). Input:
A real correlation coefficient ρ ∈ (−1, 1), a positive real boundary radius R > 0, and an
integer tracking the simulation sample size N ≥ 1000. Output: A real number representing
the estimated joint probability within the circular boundary.
 44+ Calculus Problems with Python
for Machine Learning and Quantitative Finance

1 Easy Problems
1.1 Numerical First Derivative via Forward Difference.
Computing the localized rate of change forms the core operational mechanics of checking
gradient trajectories in simple economic trends or tracking asset price velocities. Given a
continuous real-valued function f : R → R, an evaluation coordinate x, and a small step size
h, we approximate the instantaneous slope using a forward difference step. Mathematically,
this is expressed as f ′ (x) ≈ f (x+h)−f
h
(x)
. Input: A Python callable function f , a real eval-
uation point x, and a real step size h > 0. Output: A real scalar number representing the
approximated first derivative.

Examples of Different Computational Levels

def forwa rd _ d if f _ level1 (f , X , h ) : # LEVEL 1

y = []
for x in X :
df = ( f ( x + h ) - f ( x ) ) / h
y . append ( df )
return y

def forwa rd _ d if f _ level2 (f , X , h ) : # LEVEL 2

return [( f ( x + h ) - f ( x ) ) / h for x in X ]

import numpy as np

def forwa rd _ d if f _ level3 (f , X , h ) : # LEVEL 3 ( Batched Computation )

X = np . array ( X ) % Shape : (B , N ) or (N ,)
% Vectorized functional mapping over tensor arrays
return ( f ( X + h ) - f ( X ) ) / h

% Example evaluation over an array of sample nodes

X_grid = np . linspace ( -2.0 , 2.0 , 100)
df_grid = f or w a rd _ diff_level3 ( lambda x : x **2 , X_grid , 1e -5)
1.2 Numerical First Derivative via Central Difference.
Central difference formulas eliminate low-order error terms, yielding higher structural accu-
racy when approximating gradients for model regularizations or measuring instantaneous risk
sensitivities (Delta). Given a function f , an evaluation coordinate x, and a spatial step h, the
central approximation balances the interval symmetrically. Mathematically, it is defined as
f ′ (x) ≈ f (x+h)−f
2h
(x−h)
. Input: A callable function f , a real valuation point x, and a positive
real step size h > 0. Output: A real number representing the calculated derivative value.

1.3 Numerical Second Derivative Component.

The second derivative measures the local curvature or acceleration of a loss surface or pric-
ing model, serving as the fundamental metric for computing financial bond convexities or
asset Gamma. Given a function f , a target evaluation point x, and a step parameter h,
we approximate the second derivative using symmetric finite differences. Mathematically,
f ′′ (x) ≈ f (x+h)−2fh(x)+f
2
(x−h)
. Input: A callable function f , a real evaluation coordinate x,
and a small real step size h > 0. Output: A real scalar number tracking the continuous
curvature acceleration.

1.4 Riemann Sum Right-Hand Rule Quadrature.

Definite integration sums incremental changes to compute global values, such as aggregating
instantaneous continuous cash flow yields into net present value assets or calculating cu-
mulative probability distributions. Given a function f , integration limits [a, b], and a total
number of sub-intervals n, we compute
Pn the total area using right-hand boundary evaluation
anchors. Mathematically, I ≈ ∆x i=1 f (a + i∆x), where the sub-interval width is defined as
∆x = b−an . Input: A callable function f , real boundary coordinates a and b, and an integer
partition limit n ≥ 1. Output: A real scalar number representing the approximated definite
integral.

1.5 Riemann Sum Left-Hand Rule Quadrature.

Evaluating continuous accumulated areas using left-hand boundary intervals provides an al-
ternative base approximation for tracking bounding error thresholds in numerical quadrature
or expected payoff models. Given a functional equation f , integration limits [a, b], and par-
tition count
Pn−1n, we compute the total sum using left-side sample locations. Mathematically,
I ≈ ∆x i=0 f (a + i∆x) where the step size satisfies ∆x = b−a n . Input: A callable function
f , real integration parameters a and b, and an integer partition depth n ≥ 1. Output: A real
number mapping the left-hand rule area.

1.6 Average Value of a Continuous Function.

The Mean Value Theorem for Integrals implies that a continuous function achieves its long-run
average value across an interval at a specific representative node, a concept used to calculate
average volatility profiles or expected loss over uniform observation windows. Given a function
f and bounds [a, b], the average value f¯ scales the total integrated area by the inverse of the
 Rb
interval length. Mathematically, f¯ = b−a
1
a
f (x)dx. Input: A callable function f , and two
real number interval bounds a and b (a < b). Output: A real scalar number representing the
calculated continuous mean average.

1.7 Derivative of the Logistic Sigmoid Activation Function.

In deep learning, the logistic sigmoid function σ(x) = 1+e1−x maps real values to probabilities.
Evaluating its derivative is critical for updating weights during backpropagation. Mathemati-
cally, the derivative can be expressed elegantly in terms of its output: σ ′ (x) = σ(x)(1 − σ(x)).
Input: A real number x or an array of values. Output: The computed derivative value(s)
matching the input structure.

1.8 Derivative of the Rectified Linear Unit (ReLU) and LeakyReLU.

The ReLU activation function defined as f (x) = max(0, x) introduces non-linearity into ma-
chine learning networks, while LeakyReLU replaces the flat negative domain with a small
slope α to prevent dead neurons. Mathematically, the derivative of LeakyReLU is f ′ (x) = 1 if
x > 0 and f ′ (x) = α if x ≤ 0. Input: A real input coordinate x and a small leakage parameter
α (e.g., α = 0.01). Output: The evaluated derivative value.

1.9 Elasticity of a Pricing Function.

In economic modeling, algorithmic pricing, and quantitative finance, point elasticity measures
the proportional sensitivity of an output relative to an incremental change in its input variable.
Given a differentiable demand or pricing function f and an evaluation coordinate x, the
elasticity parameter scales the local first derivative by the coordinate ratio. Mathematically,
x
ϵ = f (x) · f ′ (x). Input: A callable function f , a callable first derivative function f ′ , and a
real number evaluation input x. Output: A real number representing the calculated point
elasticity score.

1.10 Linear Taylor Approximation (First-Order Tangent Matrix).

First-order Taylor expansions project linear tracking approximations near local anchor coor-
dinates, which forms the mathematical foundation for delta-hedging strategies in derivatives
trading or local linearizations of loss functions. Given a differentiable function f , a local
baseline coordinate a, and an evaluation site x, the first-order polynomial maps the local
tangent line. Mathematically, P1 (x) = f (a) + f ′ (a)(x − a). Input: A callable function f , its
derivative f ′ , an anchor point a, and an evaluation node x. Output: A real scalar estimating
the localized evaluation value.

1.11 Quadratic Taylor Expansion Polynomial.

Second-order expansions capture surface curvature to dramatically improve the fidelity of
local tracking approximations, which maps directly to delta-gamma risk models in finance
and the foundational step of Newton’s optimization method. Given a function f and its first
two derivatives f ′ , f ′′ , the local approximation includes a quadratic term centered around an
anchor coordinate a. Mathematically, P2 (x) = f (a) + f ′ (a)(x − a) + 12 f ′′ (a)(x − a)2 . Input: A
function f along with its first two derivatives f ′ and f ′′ , an anchor point a, and an evaluation
coordinate x. Output: A real scalar tracking the quadratic estimation.

1.12 Partial Derivative via Finite Difference.

Multivariate surfaces require measuring directional rates of change along individual param-
eter coordinates while holding all other input parameters constant. This forms the basis of
calculating individual weight sensitivities in machine learning. Given a multivariate func-
tion f : Rd → R, an input vector x, and a targeting dimension axis index k, we apply
a forward finite difference step along the k-th standard basis vector ek . Mathematically,
∂f f (x+hek )−f (x)
∂xk ≈ h . Input: A multivariate function f , a list of real values representing the
input vector x, an integer index k, and a real step size h > 0. Output: A real number tracking
the directional slope.

1.13 Directional Derivative Approximation.

The directional derivative projects a multi-dimensional function’s gradient along an arbitrary
directional trajectory, measuring the local rate of change along a portfolio rebalancing path or
feature-space modification vector. Given a multivariate function f , a position vector x, and
a unit direction vector v, the directional derivative is approximated by taking an incremental
step along path v. Mathematically, Dv f (x) ≈ f (x+hv)−f h
(x)
. Input: A multivariate function
f , a list representing the current position vector x, a list tracking the unit path vector v, and
a small step size h > 0. Output: A real scalar tracking the directional rate of change.

1.14 Newton-Raphson Root-Finding Iteration.

Root-finding algorithms locate the precise coordinate points where a function vanishes, a
technique used in quantitative finance to back out implied volatilities from option prices or
calculate internal rates of return (IRR). Starting from an initial guess x0 , the Newton-Raphson
algorithm updates the approximation iteratively by tracing the local tangent line down to the
x-axis. Mathematically, the update rule is defined as xk+1 = xk − ff′(x k)
(xk ) . Input: A function
′
f , its derivative f , an initial guess x0 , a real convergence tolerance ϵ, and an integer iteration
limit N . Output: A real number tracking the located root coordinate.

1.15 Binary Cross-Entropy Loss Scalar Gradient Tracker.

Binary Cross-Entropy loss is the staple objective function for classification models in machine
learning. For a single sample with true label y ∈ {0, 1} and predicted probability p = σ(z)
derived from logit z, the loss is L = −[y ln(p) + (1 − y) ln(1 − p)]. Tracking the derivative
with respect to the logit z yields the scalar error vector. Mathematically, ∂L∂z = p − y. Input:
A true label y (scalar 0 or 1) and a predicted probability scalar p ∈ (0, 1). Output: A real
scalar representing the gradient value.
2 Medium Problems
2.1 Simpson’s Rule Numerical Integration.
Simpson’s rule fits parabolic arcs across paired sub-intervals to approximate definite in-
tegrals with high accuracy, which helps evaluate continuous cumulative normal distribu-
tions used in Black-Scholes pricing models or machine learning Gaussian classifiers. Given
a function f and interval boundaries [a, b] partitioned into an even number of steps n,
the hcomposite rule weights evaluations across a non-uniform
i grid. Mathematically, I ≈
∆x
Pn/2 Pn/2−1
3 f (a) + 4 j=1 f (x2j−1 ) + 2 j=1 f (x2j ) + f (b) , where ∆x = b−a n . Input: A func-
tion f , real integration boundaries a and b, and an even integer partition count n. Output:
A real scalar number representing the integrated area.

2.2 Gradient Vector Estimation via Finite Differences.

The gradient vector aggregates all partial derivatives of a multivariate loss function into a
single directional coordinate, pointing towards the direction of maximum local ascent to guide
machine learning optimization algorithms or portfolio risk measures. Given a multivariate
function f : Rd → R and a position coordinate vector x ∈ Rd , the gradient estimator evaluates
a symmetric central difference along each coordinate axis independently. Mathematically,
h iT
∇f (x) = ∂x ∂f
1
, . . . , ∂f
∂xd
∂f
, where ∂x i
≈ f (x+hei )−f
2h
(x−hei )
. Input: A multivariate function f ,
a list representing the input vector x, and a real coordinate shift step size h > 0. Output: A
list of real numbers representing the calculated gradient vector ∇f (x).

2.3 Jacobian Matrix Construction.

The Jacobian matrix maps out the complete first-order vector derivative field for systems of
vector-valued equations, tracking how changes in multi-dimensional layer activations propa-
gate through deep neural network transformations or structural asset vector shifts. Given a
vector-valued function F : Rn → Rm , the Jacobian J ∈ Rm×n stores the partial derivatives
∂Fi
of each output component relative to each input variable: Jij = ∂x j
. We approximate these
derivatives using a forward finite difference scheme with a small coordinate shift step size
h. Input: A vector-valued function F that takes an n-dimensional vector and returns an
m-dimensional vector, a list representing the evaluation point x, and a small step size h > 0.
Output: A nested list of real numbers representing the m × n Jacobian matrix.

2.4 Hessian Matrix Estimation.

The Hessian matrix tracks the complete second-order partial derivative field of a multivari-
ate function, capturing local loss surface curvature to guide advanced optimization algo-
rithms or structure asset risk covariance representations in Risk Parity allocations. Given
a function f : Rd → R and an evaluation vector x, the Hessian matrix H ∈ Rd×d stores
2
the second-order partial derivatives: Hij = ∂x∂i ∂x
f
j
. We approximate these derivatives by
applying nested finite differences with a coordinate step size h. Mathematically, Hij ≈
f (x+hei +hej )−f (x+hei −hej )−f (x−hei +hej )+f (x−hei −hej )
4h2 .
Input: A function f , a list represent-
ing the input vector x, and a real step size h > 0. Output: A nested list of real numbers
representing the d × d Hessian matrix.

2.5 Gradient Descent Optimization Step.

Gradient descent is a foundational iterative optimization algorithm used to minimize loss
functions in machine learning and find optimal parameter configurations in quantitative trad-
ing models. Starting from an initial position vector xk , the algorithm updates the parameters
by taking a step in the direction of steepest local descent, defined by the negative gradient
−∇f (xk ), scaled by a learning rate parameter α. Mathematically, the update rule is given
by xk+1 = xk − α∇f (xk ). Input: A multivariate cost function f , an initial position vector
x0 , a learning rate scalar α > 0, an integer tracking the maximum number of iterations N ,
and a step size h for finite difference gradient estimation. Output: A list of real numbers
representing the optimized parameter vector.

2.6 Stochastic Gradient Descent (SGD) Update Loop.

Stochastic Gradient Descent accelerates optimization on massive datasets by estimating the
overall loss gradient using a single randomly selected data point or mini-batch, significantly
reducing the computational cost per iteration. P Given a global cost function decomposed into
1 M
a sum of localized component losses f (x) = M i=1 fi (x), the algorithm randomly selects an
index j ∈ {1, . . . , M } at each step and updates the parameter vector along that local negative
gradient. Mathematically, xk+1 = xk − α∇fj (xk ). Input: A list of callable loss functions
representing individual data samples, an initial parameter vector x0 , a learning rate α, and
an integer tracking the total number of update steps N . Output: A list of real numbers
representing the final optimized parameter vector.

2.7 Trapezoidal Rule with Non-Uniform Grid Spacing.

Historical financial time-series data and asset observation logs are often sampled at irregular
time intervals (tick data), requiring numerical integration methods that can handle non-
uniform grid spacing to calculate realized variance or total asset exposures. Given a sequence
of sorted domain coordinates X = [x1 , . . . , xn ] and their corresponding functional evaluations
Y = [y1 , . . . , yn ], the composite trapezoidal rule integrates the data by summing the areas
of adjacent
Pn−1 trapezoidal segments. Mathematically, the total integrated area is calculated as
I ≈ i=1 xi+12−xi (yi + yi+1 ). Input: A list of sorted real numbers X tracking coordinate
locations, and a list of real numbers Y tracking their evaluations. Output: A real scalar
number representing the total calculated integrated area.

2.8 Matrix Gradient of a Linear Regression Mean Squared Error

Objective.
In standard multivariate linear regression, the objective is to optimize the weight vector
1
w minimizing the loss function L(w) = 2N ∥Xw − y∥22 . Instead of computing derivatives
scalar by scalar, modern machine learning relies on evaluating matrix-level expressions to
update vectors simultaneously. Mathematically, the exact analytical gradient is expressed as
∇w L = N1 X T (Xw − y). Input: A NumPy design matrix X ∈ RN ×d , a target response vector
y ∈ RN , and the current weight parameter coefficients vector w ∈ Rd . Output: A NumPy
array tracking the exact analytical multi-dimensional vector gradient.

2.9 Log-Likelihood Gradient for a Univariate Gaussian Distribution.

Maximum Likelihood Estimation (MLE) underpins several foundational machine learning and
parameter calibration steps. Given an array of independent and identically distributed asset
returns or features X = [x1 , . . . , xN ], the log-likelihood function parameterized by mean µ
and variance σ 2 is optimized. Evaluating P the partial derivative with respect to µ isolates the
N
update path. Mathematically, ∂ ∂µ ln L
= σ12 i=1 (xi − µ). Input: An array of data observation
samples X, a current mean value µ, and a variance scalar σ 2 . Output: A real scalar value
specifying the exact parameter analytical gradient.

2.10 Vectorized Softmax Function Jacobian Matrix Evaluation.

The Softmax function translates an arbitrary vector of scores z ∈ Rd into a valid probability
zi
distribution where each entry si = Pe ezj . Evaluating cross-layer gradients during multi-class
categorization demands building the complete internal Jacobian matrix. Mathematically, the
entry components are defined by Jij = si (δij − sj ), where δij = 1 if i = j and 0 otherwise.
Input: A 1D array of raw classification logit scores z. Output: A 2D matrix representing the
complete calculated d × d Softmax derivative field.

2.11 Runge-Kutta 4th-Order (RK4) Differential Solver for Asset In-

terest Rates.
The 4th-order Runge-Kutta method is a highly stable iterative algorithm used to solve dif-
ferential equations, acting as a core engine for tracking continuous asset trajectories or short-
rate evolutions like the Cox-Ingersoll-Ross (CIR) interest rate framework. Given a differential
rule dydt = g(t, y) and an initial state y(t0 ) = y0 , each step evaluates four directional slopes
(k1 , k2 , k3 , k4 ) across the current interval to construct a weighted average update. Mathe-
matically, yn+1 = yn + ∆t 6 (k1 + 2k2 + 2k3 + k4 ). Input: A derivative function g(t, y), initial
coordinates (t0 , y0 ), a target endpoint tend , and a step size ∆t > 0. Output: A list of tuples
tracking the calculated solution path.

2.12 Numerical Expectation of a Continuous Payoff via Quadrature.

In automated risk management, computingR the expected payoff of a financial instrument
∞
involves evaluating the integral E[P (x)] = −∞ P (x)f (x)dx, where f (x) is the probability
density function of the underlying asset. Assuming a Gaussian distribution for f (x) truncated
across a relevant operational domain [−M, M ], we approximate the expected value using the
midpoint rule. Input: A callable payoff function P (x), mean µ and standard deviation σ
defining the density, a truncation bound M , and an integer partition count n. Output: A real
scalar representing the computed numerical expectation.
2.13 Root Finding via Secant Method for Volatility Calibration.
The secant method is an efficient iterative root-finding algorithm that avoids calculating
analytical derivatives by approximating the local slope using successive function evaluations,
making it ideal for backing out implied metrics from complex black-box derivative models.
Starting from two initial boundary guesses x0 and x1 , the algorithm iteratively updates the
root estimate by projecting the secant line down to the x-axis. Mathematically, the update
−xk−1
rule is defined as xk+1 = xk − f (xk ) f (xxkk)−f (xk−1 ) . Input: A function f , two distinct initial
guesses x0 and x1 , a real convergence tolerance ϵ, and an integer iteration limit N . Output:
A real number tracking the located root coordinate.

2.14 Richardson Extrapolation for Precision Risk Sensitivity Refine-

ment.
Richardson extrapolation is a sequence acceleration technique that combines finite difference
evaluations across multiple grid scales to systematically eliminate lower-order error terms,
significantly increasing the numerical accuracy of derivative estimates (such as checking exact
Option Gamma). Given a function f , an evaluation coordinate x, and a baseline step size h,
the algorithm combines central difference estimates computed at step sizes h and h2 . Mathe-
matically, the refined estimate is calculated using the formula D = 4·Dcentral (x,h/2)−D
3
central (x,h)
.
Input: A function f , an evaluation point x, and a baseline step size h > 0. Output: A real
scalar number representing the refined high-accuracy derivative estimate.

2.15 Total Derivative Evaluation along an Asset Price Process Tra-

jectory.
The total derivative measures the total rate of change of a multivariate value function along
a parameterized moving path, accounting for both explicit dependence on time and implicit
changes via structural underlying states. Given a portfolio valuation function f (t, S) where
asset index S(t) moves over time, the total derivative with respect to time t applies the
multivariate chain rule: df ∂f ∂f dS
dt = ∂t + ∂S dt . This can be evaluated numerically by combining
partial finite differences with spatial velocity steps. Input: A valuation function f (t, S), a
trajectory function S(t), an evaluation time t, and a small numerical step size h > 0. Output:
A real scalar number tracking the total calculated rate of change.

3 Hard Problems
3.1 Automated Backward-Mode Automatic Differentiation Engine
(Scalar Node).
Training deep neural networks and optimizing complex high-dimensional trading models re-
quires calculating exact gradients across complex computational graphs. Backward-mode
automatic differentiation computes these gradients efficiently by executing a forward pass to
calculate intermediate values, followed by a backward pass that applies the chain rule in re-
verse order to propagate derivatives from the output node back to the inputs. In this problem,
you will build a scalar computational node that maintains an explicit log of parent connections
and registers localized derivative operations, allowing it to automatically evaluate and accu-
mulate node gradients ( ∂Loss
∂x ) via a recursive .backward() call. Input: A structured sequence
of algebraic operations combining standard tracking nodes (e.g., z = ln(x) + x · y). Output:
A collection of exact scalar derivative values matching each active input node variable.

3.2 Lagrangian Multiplier Constrained Optimization Solver.

Portfolio optimization models and risk management systems frequently seek to maximize re-
turns or minimize tracking P errors subject to strict structural constraints, such as maintaining
full capital allocation ( xi = 1). The method of Lagrange multipliers converts these con-
strained problems into equivalent unconstrained systems by introducing scaling parameters λ
to penalize constraint violations. Mathematically, the objective is to locate stationary points
of the Lagrangian function L(x, λ) = f (x) + λ · g(x), which requires solving the combined sys-
tem of equations ∇x f (x) + λ∇x g(x) = 0 and g(x) = 0 using a multivariate Newton-Raphson
solver. Input: A callable objective function f (x), a constraint function g(x), an initial guess
vector [x0 , λ0 ]T , and a convergence tolerance ϵ. Output: A list of real numbers representing
the optimal constrained parameter vector x∗ .

3.3 Monte Carlo Integration over a High-Dimensional Portfolio Value-

at-Risk Hypercube.
As the number of dimensions increases, traditional grid-based numerical integration methods
suffer from the curse of dimensionality, becoming computationally intractable due to the
exponential growth of grid points. Monte Carlo integration avoids this issue by using random
sampling to approximate volumes and expectations, making it an essential tool for pricing
multi-asset options or evaluating high-dimensional risk metrics (such as multi-asset Value-
at-Risk). Given a dimension count d and a boundary radius R, the algorithm generates
N uniformly distributed points within a bounding hypercube [−R, R]d and calculates the
proportionPof points that fall within the hypervolume of the hypersphere, defined by the
d
condition i=1 x2i ≤ R2 . Input: An integer tracking total dimension size d, a real radius
parameter R > 0, and an integer specifying the total number of random sample points N .
Output: A real scalar number representing the estimated hypervolume.

3.4 Black-Scholes Partial Differential Equation (PDE) Finite Differ-

ence Solver.
The Black-Scholes PDE is a foundational continuous-time equation in financial engineering
that governs the price evolution of derivative securities based on underlying asset dynamics.
1 2 2 ∂2V
Mathematically, the option price V (S, t) satisfies the equation ∂V ∂V
∂t + 2 σ S ∂S 2 +rS ∂S −rV =
0, where σ represents volatility and r is the risk-free interest rate. This problem implements
an explicit finite difference scheme over a discrete spatial-temporal grid [0, Smax ] × [0, T ] to
backward-propagate option prices from the terminal payoff boundary V (S, T ) = max(S−K, 0)
to the initial time t = 0. Input: Execution grid parameters Smax and T , model constants
K, r, σ, and grid mesh sizes Ns and Nt . Output: A 2D array matrix containing the calculated
option prices across the asset-time grid.
3.5 Gradient Descent with Momentum Optimization Engine.
Standard gradient descent optimization paths can oscillate heavily in steep valleys or stall
at local saddle points, slowing convergence in complex machine learning landscapes. Adding
a momentum term addresses this issue by accelerating updates along consistent directions,
mimicking the physics of a marble rolling down a hill. Given an objective function f (x), the
algorithm maintains a velocity vector vk that accumulates past gradients, scaled by a decay
factor β, and updates the parameter vector along this momentum-adjusted path. Mathemat-
ically, the update rules are defined as vk+1 = βvk + α∇f (xk ) and xk+1 = xk − vk+1 . Input: A
cost function f , an initial parameter vector x0 , a learning rate α, a momentum decay scalar
β ∈ [0, 1), and a maximum iteration cap N . Output: A list of real numbers representing the
optimized parameter vector.

3.6 Backpropagation Gradient through a Single Hidden Layer Neu-

ral Network.
To train artificial networks, errors must be analytically propagated backward from the final
objective layer through matrix transformations. Given an input matrix X, weights W1 and
W2 , and a bias vector b1 , the network predicts ŷ = σ(W2 · ReLU(W1 · X + b1 )). Your task is to
calculate the explicit analytical multivariate gradient matrix ∂Loss
∂W1 for a Mean Squared Error
loss layer using matrix calculus chains and element-wise masking operations. Input: Input
data matrix X, initial transformation weight matrices W1 and W2 , target matrix labels y,
and bias parameters. Output: A multi-dimensional array matrix tracing the calculated layer
parameter updates.

3.7 Secant-Based Quasi-Newton BFGS Optimization Engine Step.

The Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm is a powerful quasi-Newton opti-
mization method that avoids the high computational cost of calculating the full Hessian matrix
at each step by iteratively updating an approximation of the inverse Hessian matrix Hk . Given
a displacement vector sk = xk+1 − xk and a change in gradient yk = ∇f (xk+1 ) − ∇f (xk ), the
inverse Hessian matrix is updated using a rank-two correction formula that preserves symme-
try and positive definiteness. Mathematically, Hk+1 = (I − ρk sk ykT )Hk (I − ρk yk sTk ) + ρk sk sTk ,
where ρk = yT1sk . Input: A current position vector xk , a gradient vector ∇f (xk ), the current
k
inverse Hessian approximation matrix Hk , and a line search step size α. Output: A tuple
containing the updated position vector xk+1 and the updated inverse Hessian matrix Hk+1 .

3.8 Numerical Inversion of a Characteristic Function via Integra-

tion.
Advanced asset valuation models (such as the Heston stochastic volatility model) define op-
tion pricing equations in the frequency domain through characteristic functions, requiring
numerical inversion via Fourier transform
 −iϕ ln(K) integrals.
 The contract value involves evaluating
R∞
integrals of the form I = 0 Re e iϕ
ψ(ϕ)
dϕ. This problem implements an automated
high-resolution composite adaptive integration engine to solve this pricing term over a trun-
cated calculation window. Input: A complex-valued callable characteristic function ψ, a
target strike threshold K, and an integration limit M . Output: A real scalar tracking the
exact inverted valuation probability.

3.9 Ornstein-Uhlenbeck Mean-Reverting Process Drift and Diffu-

sion Numerical Matcher.
The Ornstein-Uhlenbeck (OU) process is defined by the stochastic differential equation dxt =
θ(µ − xt )dt + σdWt , forming the foundation for pairs trading and volatility modeling. Cali-
brating the system parameters from discrete time-series historical observations relies on maxi-
mizing the continuous transition log-likelihood. You will construct a solver that takes discrete
observations and tracks the numerical calculus gradient optimization path to map the param-
eters (θ, µ, σ 2 ). Input: An array of historically tracked asset values X and observation time
interval ∆t. Output: A vector of optimized model coefficients matching drift and diffusion
profiles.

3.10 Value-at-Risk (VaR) and Expected Shortfall Continuous Gra-

dient Tracking.
Quantifying portfolio risk requires monitoring continuous structural performance indices like
Expected Shortfall (ES). When individual asset return densities are modeled as smooth joint
mixtures, the boundary point defining Value-at-Risk (VaR) moves continuously based on
portfolio allocations w. Your engine will compute the numerical multivariate gradient of
the expected shortfall boundary area with respect to the weight allocations w to enable
downstream risk-managed portfolio optimization loops. Input: An allocation vector w, a
multivariate asset returns dataset R, and a target tail risk confidence threshold α (e.g., α =
0.95). Output: A risk sensitivity vector tracking allocation updates.

3.11 Implicit Crank-Nicolson Scheme for Finite Difference Asset

Diffusion.
The Crank-Nicolson method is a numerically stable implicit finite difference scheme used to
solve partial differential equations like option pricing or heat diffusion frameworks, avoiding
∂2u
the strict step-size limitations of explicit methods. Given a diffusion rule ∂u
∂t = κ ∂x2 , the algo-
rithm forms an implicit matrix system at each step by averaging a forward step in time with
a central step in space. This updates the state across the spatial grid by solving a tridiagonal
linear equation at each chronological step. Input: Initial spatial grid profile array U0 , diffu-
sion coefficient constant κ, boundary condition constants, and grid parameters ∆x, ∆t, Nt .
Output: A 2D array matrix tracking the complete evolution of the spatial-temporal diffusion
profile.

3.12 Continuous Gradient Flow Trajectory Tracker for Deep Loss

Optimization.
Gradient flow models the continuous-time limit of gradient descent, tracing an idealized
smooth parameter path that moves along the direction of steepest local descent to provide
theoretical insights into optimization dynamics. Mathematically, the parameter trajectory
satisfies the system of ordinary differential equations dx dt = −∇f (x(t)), with an initial con-
dition x(0) = x0 . We track this continuous optimization path numerically by solving the
system of equations using a high-order adaptive differential solver. Input: A multivariate
cost function f, an initial position vector x0 , a total time horizon T , and an integration time
step size dt. Output: A nested list tracking the sequential coordinate positions along the
continuous gradient flow trajectory.

3.13 Adam (Adaptive Moment Estimation) Optimizer Core Engine

Step.
The Adam optimization algorithm is the industry standard for minimizing complex empirical
loss functions when training deep neural networks. It maintains running tracking averages
of both the first moment (gradients) and the second raw moment (uncentered squared gra-
dients) to dynamically adjust learning rates parameter-by-parameter. Mathematically, the
mt vt
system incorporates bias-correction transformations: m̂t = 1−β t and v̂t = 1−β t , yielding
1 2
α
parameter updates xt+1 = xt − √v̂ +ϵ m̂t . Input: Parameter vector x, local gradient vector
t
g, historical moment tracking structures, hyperparameters (α, β1 , β2 , ϵ), and the integer step
count t. Output: A tuple containing the updated parameter vector and refreshed moment
tracking structures.

3.14 Bivariate Newton-Raphson Optimization Engine.

Locating the precise stationary points of a multi-variable surface requires finding the coordi-
nate locations where all partial derivatives vanish simultaneously, a key step for solving multi-
asset equilibrium conditions or maximizing joint likelihood functions. Given a two-variable
function f (x, y), the bivariate Newton-Raphson algorithm updates the coordinate position
vector by adjusting it along the direction of the inverse Hessian matrix scaled by the local
gradient vector. Mathematically, the update rule is defined as Xk+1 = Xk −H −1 (Xk )∇f (Xk ),
which we evaluate numerically by constructing and inverting the local second-order finite dif-
ference matrices at each step. Input: A two-variable objective function f (x, y), an initial
position vector [x0 , y0 ]T , a convergence tolerance threshold ϵ, and an iteration limit N . Out-
put: A list of real numbers tracking the located stationary coordinate vector.

3.15 Numerical Inversion of the Cumulative Distribution Function

for Copula Modeling.
In structural asset correlation and multi-variable risk management, financial copulas isolate
dependencies by projecting marginal features into standard uniform spaces. This process
relies on the numerical inversion of a distribution’s CDF, solving F (x) = u for x given a
probability level u ∈ (0, 1). You will build an engine that evaluates the forward continu-
ous CDF via numerical quadrature steps and applies an internal Newton-Raphson boundary
search iteration to locate high-precision inverse quantile points. Input: A callable continuous
probability density function (PDF) f , a target uniform probability density coordinate u, and
an initial root estimation anchor. Output: A real scalar tracking the exact computed quantile
value.