সিঙ্গুলার ভ্যালু ডিকম্পোজিশন (SVD)

import math

# A হলো 3x2 ম্যাট্রিক্স (৩টি সারি, ২টি কলাম) -- একটি ইলাস্ট্রেটিভ উদাহরণ
A = [
    [3.0, 1.0],
    [1.0, 3.0],
    [1.0, 1.0],
]

def transpose(M):
    rows, cols = len(M), len(M[0])
    return [[M[r][c] for r in range(rows)] for c in range(cols)]

def matmul(M1, M2):
    r1, c1 = len(M1), len(M1[0])
    r2, c2 = len(M2), len(M2[0])
    result = [[0.0]*c2 for _ in range(r1)]
    for i in range(r1):
        for j in range(c2):
            s = 0.0
            for k in range(c1):
                s += M1[i][k]*M2[k][j]
            result[i][j] = s
    return result

def matvec(M, v):
    return [sum(M[i][j]*v[j] for j in range(len(v))) for i in range(len(M))]

def vec_norm(v):
    return math.sqrt(sum(x*x for x in v))

def normalize(v):
    n = vec_norm(v)
    return [x/n for x in v]

At = transpose(A)
AtA = matmul(At, A)
print("A^T A =", AtA)
print()

# M8-এর পাওয়ার মেথড, এখানে A^T A-এর উপর প্রয়োগ করা হচ্ছে
v = [1.0, 0.0]
print(f"{'ইটারেশন':>8} | {'v (normalized)':<26} | {'lambda estimate':>16}")
lambda_est = None
for i in range(1, 13):
    w = matvec(AtA, v)
    lambda_est = vec_norm(w)
    v = normalize(w)
    print(f"{i:8d} | {str([round(x,6) for x in v]):<26} | {lambda_est:16.8f}")

lambda1 = lambda_est
sigma1 = math.sqrt(lambda1)
v1 = v
print()
print(f"কনভার্জড ডমিনেন্ট আইগেনভ্যালু (lambda1): {lambda1:.6f}")
print(f"রাইট সিঙ্গুলার ভেক্টর v1: {[round(x,6) for x in v1]}")
print(f"প্রথম সিঙ্গুলার ভ্যালু sigma1 = sqrt(lambda1): {sigma1:.6f}")

# u1 = (1/sigma1) * A v1
Av1 = matvec(A, v1)
u1 = [x/sigma1 for x in Av1]
print(f"লেফট সিঙ্গুলার ভেক্টর u1 = (1/sigma1) A v1: {[round(x,6) for x in u1]}")
print(f"u1-এর নর্ম (~1 হওয়া উচিত): {vec_norm(u1):.6f}")

# rank-1 রিকনস্ট্রাকশন sigma1 * u1 * v1^T
recon = [[sigma1*u1[i]*v1[j] for j in range(2)] for i in range(3)]
print()
print("Rank-1 রিকনস্ট্রাকশন sigma1*u1*v1^T:")
for row in recon:
    print([round(x, 6) for x in row])
print("মূল A:")
for row in A:
    print(row)

diff = [[A[i][j]-recon[i][j] for j in range(2)] for i in range(3)]
frob_err = math.sqrt(sum(x*x for row in diff for x in row))
print()
print(f"রেসিডুয়ালের Frobenius নর্ম (A - rank-1 reconstruction): {frob_err:.6f}")