z
"""
Resources:
- https://en.wikipedia.org/wiki/Definite_symmetric_matrix
"""
from typing import Any

import numpy as np

def _is_matrix_spd(matrix: np.ndarray) -> bool:
"""
Returns True if input matrix is symmetric positive definite.
Returns False otherwise.

For a matrix to be SPD, all eigenvalues must be positive.

>>> import numpy as np
>>> matrix = np.array([
... [4.12401784, -5.01453636, -0.63865857],
... [-5.01453636, 12.33347422, -3.40493586],
... [-0.63865857, -3.40493586,  5.78591885]])
>>> _is_matrix_spd(matrix)
True
>>> matrix = np.array([
... [0.34634879,  1.96165514,  2.18277744],
... [0.74074469, -1.19648894, -1.34223498],
... [-0.7687067 ,  0.06018373, -1.16315631]])
>>> _is_matrix_spd(matrix)
False
"""
# Ensure matrix is square.
assert np.shape(matrix)[0] == np.shape(matrix)[1]

# If matrix not symmetric, exit right away.
if np.allclose(matrix, matrix.T) is False:
return False

# Get eigenvalues and eignevectors for a symmetric matrix.
eigen_values, _ = np.linalg.eigh(matrix)

# Check sign of all eigenvalues.
# np.all returns a value of type np.bool_
return bool(np.all(eigen_values > 0))

def _create_spd_matrix(dimension: int) -> Any:
"""
Returns a symmetric positive definite matrix given a dimension.

Input:
dimension gives the square matrix dimension.

Output:
spd_matrix is an diminesion x dimensions symmetric positive definite (SPD) matrix.

>>> import numpy as np
>>> dimension = 3
>>> spd_matrix = _create_spd_matrix(dimension)
>>> _is_matrix_spd(spd_matrix)
True
"""
random_matrix = np.random.randn(dimension, dimension)
spd_matrix = np.dot(random_matrix, random_matrix.T)
assert _is_matrix_spd(spd_matrix)
return spd_matrix

spd_matrix: np.ndarray,
max_iterations: int = 1000,
tol: float = 1e-8,
) -> Any:
"""
Returns solution to the linear system np.dot(spd_matrix, x) = b.

Input:
spd_matrix is an NxN Symmetric Positive Definite (SPD) matrix.

Output:
x is an Nx1 vector that is the solution vector.

>>> import numpy as np
>>> spd_matrix = np.array([
... [8.73256573, -5.02034289, -2.68709226],
... [-5.02034289,  3.78188322,  0.91980451],
... [-2.68709226,  0.91980451,  1.94746467]])
>>> b = np.array([
... [-5.80872761],
... [ 3.23807431],
... [ 1.95381422]])
array([[-0.63114139],
[-0.01561498],
[ 0.13979294]])
"""
# Ensure proper dimensionality.
assert np.shape(spd_matrix)[0] == np.shape(spd_matrix)[1]
assert _is_matrix_spd(spd_matrix)

# Initialize solution guess, residual, search direction.
p0 = np.copy(r0)

# Set initial errors in solution guess and residual.
error_residual = 1e9
error_x_solution = 1e9
error = 1e9

# Set iteration counter to threshold number of iterations.
iterations = 0

while error > tol:
# Save this value so we only calculate the matrix-vector product once.
w = np.dot(spd_matrix, p0)

# The main algorithm.

# Update search direction magnitude.
alpha = np.dot(r0.T, r0) / np.dot(p0.T, w)
# Update solution guess.
x = x0 + alpha * p0
# Calculate new residual.
r = r0 - alpha * w
# Calculate new Krylov subspace scale.
beta = np.dot(r.T, r) / np.dot(r0.T, r0)
# Calculate new A conjuage search direction.
p = r + beta * p0

# Calculate errors.
error_residual = np.linalg.norm(r - r0)
error_x_solution = np.linalg.norm(x - x0)
error = np.maximum(error_residual, error_x_solution)

# Update variables.
x0 = np.copy(x)
r0 = np.copy(r)
p0 = np.copy(p)

# Update number of iterations.
iterations += 1
if iterations > max_iterations:
break

return x

"""
>>> test_conjugate_gradient()  # self running tests
"""
# Create linear system with SPD matrix and known solution x_true.
dimension = 3
spd_matrix = _create_spd_matrix(dimension)
x_true = np.random.randn(dimension, 1)
b = np.dot(spd_matrix, x_true)

# Numpy solution.
x_numpy = np.linalg.solve(spd_matrix, b)

# Our implementation.