import time
import numpy as np
from math import sqrt, cos, acos, pi
A = np.array([
[0, 1, 1],
[sqrt(2), 2, 0],
[0, 1, 1]
])
def compute_trace_manual(A):
trace = 0
for i in range(len(A)):
trace += A[i][i]
return trace
def compute_determinant_manual(A):
a = A[0][0]; b = A[0][1]; c = A[0][2]
d = A[1][0]; e = A[1][1]; f = A[1][2]
g = A[2][0]; h = A[2][1]; i = A[2][2]
determinant = (a * (e * i - f * h) -
b * (d * i - f * g) +
c * (d * h - e * g))
return determinant
def compute_characteristic_polynomial_coefficients(A):
p = compute_trace_manual(A)
minor1 = A[1][1]*A[2][2] - A[1][2]*A[2][1]
minor2 = A[0][0]*A[2][2] - A[0][2]*A[2][0]
minor3 = A[0][0]*A[1][1] - A[0][1]*A[1][0]
q = minor1 + minor2 + minor3
r = compute_determinant_manual(A)
return [1, -p, q, -r]
def solve_cubic(coeffs):
a, b, c, d = coeffs
p = (3*a*c - b**2) / (3*a**2)
q = (2*b**3 - 9*a*b*c + 27*a**2*d) / (27*a**3)
discriminant = (q/2)**2 + (p/3)**3
if discriminant > 0:
u = (-q/2 + sqrt(discriminant))**(1/3)
v = (-q/2 - sqrt(discriminant))**(1/3)
root1 = u + v - b / (3*a)
roots = [root1]
elif discriminant == 0:
u = (-q/2)**(1/3)
root1 = 2*u - b / (3*a)
root2 = -u - b / (3*a)
roots = [root1, root2]
else:
r = sqrt(-p**3/27)
theta = acos(-q/(2*r))
r = (-p/3)**(1/2)
root1 = 2*r*cos(theta/3) - b / (3*a)
root2 = 2*r*cos((theta + 2*pi)/3) - b / (3*a)
root3 = 2*r*cos((theta + 4*pi)/3) - b / (3*a)
roots = [root1, root2, root3]
return roots
def compute_eigenvalues_eigenvectors_manual(A):
coeffs = compute_characteristic_polynomial_coefficients(A)
eigenvalues = solve_cubic(coeffs)
eigenvectors = []
for eigenvalue in eigenvalues:
M = A - eigenvalue * np.identity(3)
row1 = M[0]
row2 = M[1]
eigenvector = np.cross(row1, row2)
if np.linalg.norm(eigenvector) == 0:
row1 = M[1]
row2 = M[2]
eigenvector = np.cross(row1, row2)
if np.linalg.norm(eigenvector) == 0:
row1 = M[0]
row2 = M[2]
eigenvector = np.cross(row1, row2)
eigenvector = eigenvector / np.linalg.norm(eigenvector)
eigenvectors.append(eigenvector.real)
eigenvalues = np.array(eigenvalues).real
eigenvectors = np.array(eigenvectors).T
return eigenvalues, eigenvectors
start_time = time.time()
custom_eigenvalues, custom_eigenvectors = compute_eigenvalues_eigenvectors_manual(A)
custom_runtime = time.time() - start_time
start_time = time.time()
numpy_eigenvalues, numpy_eigenvectors = np.linalg.eig(A)
numpy_runtime = time.time() - start_time
print("Custom Eigenvalues:\n", custom_eigenvalues)
print("Custom Eigenvectors:\n", custom_eigenvectors)
print("Custom Function Runtime:", custom_runtime, "seconds\n")
print("NumPy Eigenvalues:\n", numpy_eigenvalues)
print("NumPy Eigenvectors:\n", numpy_eigenvectors)
print("NumPy Function Runtime:", numpy_runtime, "seconds")
Custom Eigenvalues: [ 2.79004402e+00 -6.66133815e-16 2.09955984e-01] Custom Eigenvectors: [[ 0.43834959 -0.70710678 -0.61731105] [ 0.78466507 0.5 0.4877029 ] [ 0.43834959 -0.5 -0.61731105]] Custom Function Runtime: 0.003560781478881836 seconds NumPy Eigenvalues: [ 2.79004402e+00 -2.63339565e-19 2.09955984e-01] NumPy Eigenvectors: [[-0.43834959 0.70710678 0.61731105] [-0.78466507 -0.5 -0.4877029 ] [-0.43834959 0.5 0.61731105]] NumPy Function Runtime: 0.0005056858062744141 seconds
def compute_covariance_matrix_manual(A):
n_samples = A.shape[0]
mean_centered = A - np.mean(A, axis=0)
covariance_matrix = np.zeros((A.shape[1], A.shape[1]))
for i in range(A.shape[1]):
for j in range(A.shape[1]):
covariance_matrix[i][j] = np.sum(mean_centered[:, i] * mean_centered[:, j]) / (n_samples - 1)
return covariance_matrix
def compute_pca_manual(A, n_components=None):
A_meaned = A - np.mean(A, axis=0)
covariance_matrix = compute_covariance_matrix_manual(A)
#Here, We've reused our manually implemented eigenvalue and eigenvector computation function.
eigenvalues, eigenvectors = compute_eigenvalues_eigenvectors_manual(covariance_matrix)
sorted_idx = np.argsort(eigenvalues)[::-1]
sorted_eigenvalues = eigenvalues[sorted_idx]
sorted_eigenvectors = eigenvectors[:, sorted_idx]
explained_variance = sorted_eigenvalues / np.sum(sorted_eigenvalues)
if n_components is not None:
sorted_eigenvectors = sorted_eigenvectors[:, :n_components]
explained_variance = explained_variance[:n_components]
principal_components = np.dot(A_meaned, sorted_eigenvectors)
return principal_components, explained_variance, sorted_eigenvectors
start_time = time.time()
custom_pcs, custom_explained_variance, custom_eigenvectors = compute_pca_manual(A)
custom_pca_runtime = time.time() - start_time
from sklearn.decomposition import PCA
start_time = time.time()
pca = PCA()
pca.fit(A)
sklearn_pca_runtime = time.time() - start_time
print("\nCustom Principal Components:\n", custom_pcs)
print("Explained Variance:\n", custom_explained_variance)
print("Custom PCA Runtime:", custom_pca_runtime, "seconds\n")
print("Scikit-Learn Principal Components:\n", pca.transform(A))
print("Explained Variance:\n", pca.explained_variance_ratio_)
print("Scikit-Learn PCA Runtime:", sklearn_pca_runtime, "seconds")
Custom Principal Components: [[ 0.66666667] [-1.33333333] [ 0.66666667]] Explained Variance: [1.] Custom PCA Runtime: 0.0015196800231933594 seconds Scikit-Learn Principal Components: [[-6.66666667e-01 0.00000000e+00 0.00000000e+00] [ 1.33333333e+00 2.22044605e-16 1.11022302e-16] [-6.66666667e-01 0.00000000e+00 0.00000000e+00]] Explained Variance: [1.00000000e+00 1.42738564e-32 7.48397169e-34] Scikit-Learn PCA Runtime: 0.0016303062438964844 seconds
def matrix_multiply_manual(X, Y):
result = np.zeros((X.shape[0], Y.shape[1]))
for i in range(X.shape[0]):
for j in range(Y.shape[1]):
sum = 0
for k in range(X.shape[1]):
sum += X[i][k] * Y[k][j]
result[i][j] = sum
return result
def compute_svd_manual(A):
At = A.T
AtA = matrix_multiply_manual(At, A)
AAt = matrix_multiply_manual(A, At)
#Here, We've reused our manually implemented eigenvalue and eigenvector computation function.
eigenvalues_V, eigenvectors_V = compute_eigenvalues_eigenvectors_manual(AtA)
eigenvalues_U, eigenvectors_U = compute_eigenvalues_eigenvectors_manual(AAt)
singular_values = np.sqrt(eigenvalues_V)
sorted_idx = np.argsort(singular_values)[::-1]
singular_values = singular_values[sorted_idx]
U = eigenvectors_U[:, sorted_idx]
V = eigenvectors_V[:, sorted_idx]
Sigma = np.zeros_like(A, dtype=float)
min_dim = min(A.shape)
for i in range(min_dim):
Sigma[i][i] = singular_values[i]
return U, Sigma, V.T
start_time = time.time()
custom_U, custom_Sigma, custom_Vt = compute_svd_manual(A)
custom_svd_runtime = time.time() - start_time
start_time = time.time()
numpy_U, numpy_S, numpy_Vt = np.linalg.svd(A)
numpy_svd_runtime = time.time() - start_time
print("\nCustom SVD U:\n", custom_U)
print("Custom Singular Values:\n", np.diag(custom_Sigma))
print("Custom SVD V^T:\n", custom_Vt)
print("Custom SVD Runtime:", custom_svd_runtime, "seconds\n")
print("NumPy SVD U:\n", numpy_U)
print("NumPy Singular Values:\n", numpy_S)
print("NumPy SVD V^T:\n", numpy_Vt)
print("NumPy SVD Runtime:", numpy_svd_runtime, "seconds")
Custom SVD U: [[ 0.40824829 -0.57735027 -0.70710678] [ 0.81649658 0.57735027 0. ] [ 0.40824829 -0.57735027 0.70710678]] Custom Singular Values: [2.82842712 1.41421356 0. ] Custom SVD V^T: [[ 4.08248290e-01 8.66025404e-01 2.88675135e-01] [ 5.77350269e-01 -3.62597321e-16 -8.16496581e-01] [ 7.07106781e-01 -5.00000000e-01 5.00000000e-01]] Custom SVD Runtime: 0.004405021667480469 seconds NumPy SVD U: [[-0.40824829 0.57735027 -0.70710678] [-0.81649658 -0.57735027 0. ] [-0.40824829 0.57735027 0.70710678]] NumPy Singular Values: [2.82842712 1.41421356 0. ] NumPy SVD V^T: [[-4.08248290e-01 -8.66025404e-01 -2.88675135e-01] [-5.77350269e-01 -4.89172797e-17 8.16496581e-01] [ 7.07106781e-01 -5.00000000e-01 5.00000000e-01]] NumPy SVD Runtime: 0.0004723072052001953 seconds