স্ট্রাসেনের গুণন ও মিডিয়ান ফাইন্ডিং অ্যালগরিদম

১ · সমস্যা — বড় ম্যাট্রিক্স গুণন

দুটি $n \times n$ ম্যাট্রিক্স $A$ ও $B$ গুণ করে $C = AB$ বের করার নেইভ পদ্ধতি হলো প্রতিটি $c_{ij} = \sum_{k} a_{ik} b_{kj}$ সরাসরি হিসাব করা — তিনটি নেস্টেড লুপ, $\Theta(n^3)$ সময়। প্রশ্ন হলো: D&C দিয়ে কি এর চেয়ে ভালো করা সম্ভব?

সরল ব্লক-ভিত্তিক D&C পদ্ধতি চেষ্টা করলে: $A$ ও $B$-কে চারটি $n/2 \times n/2$ ব্লকে ভাগ করে ($A_{11}, A_{12}, A_{21}, A_{22}$ এবং $B_{11}, B_{12}, B_{21}, B_{22}$), সাধারণ ম্যাট্রিক্স গুণনের নিয়মে:

এখানে মোট $8$টি ব্লক-গুণন লাগে (প্রতিটি $C_{ij}$-এ ২টি করে) এবং কিছু $O(n^2)$ যোগ। রিকারেন্স $T(n) = 8T(n/2) + O(n^2)$ — কিন্তু $n^{\log_2 8} = n^3$, তাই মাস্টার থিওরেমের কেস ১ প্রয়োগ করলে $T(n) = \Theta(n^3)$ — নেইভ পদ্ধতির চেয়ে কোনো লাভ নেই! এখানেই স্ট্রাসেনের অন্তর্দৃষ্টি কাজে আসে।

২ · স্ট্রাসেনের মূল অন্তর্দৃষ্টি — ৮টির বদলে ৭টি গুণন

ম্যাট্রিক্স গুণন ($O(n^3)$-ঘেঁষা, ব্যয়বহুল) আর যোগ/বিয়োগ ($O(n^2)$, সস্তা) — এই দুইয়ের খরচে বিশাল ফারাক। স্ট্রাসেন দেখান যে সাতটি চতুরভাবে নির্মিত গুণফল দিয়েই (কিছু অতিরিক্ত যোগ-বিয়োগের বিনিময়ে) সবগুলো $C_{ij}$ ব্লক পুনর্গঠন করা সম্ভব — অষ্টম গুণনটি সম্পূর্ণ বাদ দেওয়া যায়:

এবং চূড়ান্ত ব্লকগুলো:

এই সূত্রগুলো সরাসরি বীজগণিতের মাধ্যমে যাচাই করা যায় — উদাহরণস্বরূপ, $M_1 + M_4 - M_5 + M_7$ খুলে সমস্ত পদ বিস্তার করলে দেখা যায় প্রায় সব ক্রস-টার্ম একে অপরকে কেটে দেয়, শুধু $A_{11}B_{11} + A_{12}B_{21}$ অবশিষ্ট থাকে — যা ঠিক $C_{11}$-এর সংজ্ঞা। নিচের কোড সেলে এই সূত্রগুলো হুবহু ইমপ্লিমেন্ট করে নেইভ গুণনের সাথে মিলিয়ে এই দাবিটি সরাসরি যাচাই করা হচ্ছে — শুধু বিশ্বাস করে নয়।

৩ · রিকারেন্স ডেরিভেশন — $\Theta(n^{\log_2 7})$

প্রতিটি $M_i$ একটি $n/2 \times n/2$ ম্যাট্রিক্স গুণন (রিকার্সিভ কল), এবং যোগ/বিয়োগ/কমবাইন মিলিয়ে $O(n^2)$ অতিরিক্ত কাজ। তাই রিকারেন্স:

মাস্টার থিওরেম প্রয়োগ করতে ($a=7$, $b=2$, $f(n) = \Theta(n^2)$): প্রথমে $n^{\log_b a}$ হিসাব করা হয়:

যেহেতু $f(n) = \Theta(n^2)$ এবং $2 < \log_2 7 \approx 2.807$, তাই $f(n) = O(n^{\log_2 7 - \epsilon})$ কোনো $\epsilon > 0$-এর জন্য (এখানে $\epsilon \approx 0.807$) — এটি মাস্টার থিওরেমের কেস ১ (L13 দ্রষ্টব্য)। কেস ১ বলে:

অর্থাৎ স্ট্রাসেনের অ্যালগরিদম $\Theta(n^{2.807})$ সময়ে চলে — নেইভ $\Theta(n^3)$-এর চেয়ে অ্যাসিম্পটোটিকভাবে দ্রুত (যদিও ব্যবহারিকভাবে ছোট ম্যাট্রিক্সে কনস্ট্যান্ট ফ্যাক্টর ও অতিরিক্ত মেমরি-ব্যবহারের কারণে নেইভ পদ্ধতিই দ্রুত হতে পারে — এটি একটি সাধারণ প্যাটার্ন যা M13-এ ("অ্যালগরিদম অ্যানালাইসিসের সাধারণ ভুল") আবার আসবে)।

৪ · মিডিয়ান ফাইন্ডিং — নেইভ বনাম worst-case $O(n)$ সিলেক্ট

একটি সম্পূর্ণ ভিন্ন সমস্যা, কিন্তু একই D&C পরিবারের: একটি আনসর্টেড অ্যারের $k$-তম ক্ষুদ্রতম এলিমেন্ট (যেমন মিডিয়ান) বের করা। নেইভ পদ্ধতি — পুরো অ্যারে সর্ট করে $k$-তম ইনডেক্সে দেখা — $\Theta(n \log n)$ সময় নেয়, যদিও আমরা পুরো সর্টেড ক্রম চাই না, শুধু একটি এলিমেন্ট।

মিডিয়ান-অফ-মিডিয়ানস (Select) অ্যালগরিদম worst-case $O(n)$-এ এটি অর্জন করে একটি চতুর পিভট-বাছাই কৌশল দিয়ে: অ্যারেকে ৫টি করে গ্রুপে ভাগ করা হয়, প্রতিটি গ্রুপের মিডিয়ান বের করা হয় (ছোট, তাই $O(1)$ প্রতি গ্রুপ), তারপর সেই মিডিয়ানগুলোরও মিডিয়ান রিকার্সিভভাবে বের করে সেটিকে পিভট হিসেবে ব্যবহার করা হয়। এই বিশেষ পিভট নিশ্চিত করে যে পার্টিশনের প্রতিটি ভাগে অন্তত প্রায় $3n/10$টি এলিমেন্ট থাকবে (কখনোই খুব বেশি অসামঞ্জস্যপূর্ণ স্প্লিট হবে না), যা রিকারেন্স $T(n) \le T(n/5) + T(7n/10) + O(n)$ দেয়। যেহেতু $\frac{1}{5} + \frac{7}{10} = \frac{9}{10} < 1$ (উপ-প্রবলেমের আকারের যোগফল সবসময় মূল আকারের একটি ধ্রুবক ভগ্নাংশের কম থাকে), এই রিকারেন্স $T(n) = O(n)$-এ সমাধান হয় (একটি জ্যামিতিক ধারার মতো ক্ষয়িষ্ণু সিরিজ যোগ করে — সম্পূর্ণ প্রমাণ এখানে সংক্ষিপ্ত রাখা হলো, কারণ এই কোর্সের মূল লক্ষ্য D&C টেমপ্লেট চেনা, প্রতিটি নির্বাচন-অ্যালগরিদমের প্রতিটি ধ্রুবক নিখুঁতভাবে অপ্টিমাইজ করা নয়)।

৫ · কম্পিউটেশনাল যাচাই — স্ট্রাসেনের গুণন

নিচের কোডে স্ট্রাসেনের অ্যালগরিদম $4 \times 4$ (এবং $2 \times 2$) ম্যাট্রিক্সের জন্য ইমপ্লিমেন্ট করা হয়েছে — উপরের সূত্রগুলো হুবহু কোডে রূপান্তরিত — এবং এর আউটপুট একটি সরল নেইভ ট্রিপল-লুপ গুণনের সাথে একাধিক র‍্যান্ডম ম্যাট্রিক্সে নিখুঁতভাবে মিলিয়ে দেখা হয়েছে।

import random

def matrix_add(A, B):
    return [[A[i][j] + B[i][j] for j in range(len(A[0]))] for i in range(len(A))]

def matrix_sub(A, B):
    return [[A[i][j] - B[i][j] for j in range(len(A[0]))] for i in range(len(A))]

def split_quadrants(M):
    n = len(M)
    mid = n // 2
    A11 = [row[:mid] for row in M[:mid]]
    A12 = [row[mid:] for row in M[:mid]]
    A21 = [row[:mid] for row in M[mid:]]
    A22 = [row[mid:] for row in M[mid:]]
    return A11, A12, A21, A22

def combine_quadrants(C11, C12, C21, C22):
    top = [r1 + r2 for r1, r2 in zip(C11, C12)]
    bottom = [r1 + r2 for r1, r2 in zip(C21, C22)]
    return top + bottom

def strassen_multiply(A, B):
    n = len(A)
    if n == 1:
        return [[A[0][0] * B[0][0]]]

    A11, A12, A21, A22 = split_quadrants(A)
    B11, B12, B21, B22 = split_quadrants(B)

    M1 = strassen_multiply(matrix_add(A11, A22), matrix_add(B11, B22))
    M2 = strassen_multiply(matrix_add(A21, A22), B11)
    M3 = strassen_multiply(A11, matrix_sub(B12, B22))
    M4 = strassen_multiply(A22, matrix_sub(B21, B11))
    M5 = strassen_multiply(matrix_add(A11, A12), B22)
    M6 = strassen_multiply(matrix_sub(A21, A11), matrix_add(B11, B12))
    M7 = strassen_multiply(matrix_sub(A12, A22), matrix_add(B21, B22))

    C11 = matrix_add(matrix_sub(matrix_add(M1, M4), M5), M7)
    C12 = matrix_add(M3, M5)
    C21 = matrix_add(M2, M4)
    C22 = matrix_add(matrix_sub(matrix_add(M1, M3), M2), M6)

    return combine_quadrants(C11, C12, C21, C22)

def naive_multiply(A, B):
    n, m, k = len(A), len(B[0]), len(B)
    C = [[0] * m for _ in range(n)]
    for i in range(n):
        for j in range(m):
            total = 0
            for x in range(k):
                total += A[i][x] * B[x][j]
            C[i][j] = total
    return C

def random_matrix(rng, n):
    return [[rng.randint(-9, 9) for _ in range(n)] for _ in range(n)]

all_match = True
for size in [2, 4]:
    for seed in [1, 2, 3, 4, 5]:
        rng = random.Random(seed)
        A = random_matrix(rng, size)
        B = random_matrix(rng, size)
        result_strassen = strassen_multiply(A, B)
        result_naive = naive_multiply(A, B)
        match = (result_strassen == result_naive)
        all_match = all_match and match
        print(f"size={size} seed={seed}: স্ট্রাসেন == নেইভ ? {match}")

print(f"\nসব কেসে স্ট্রাসেন ও নেইভ গুণন মিলেছে: {all_match}")
assert all_match

৬ · কম্পিউটেশনাল যাচাই — Quickselect দিয়ে মিডিয়ান

import random

def quickselect(a, k, rng):
    """a-এর মধ্যে k-তম ক্ষুদ্রতম এলিমেন্ট (0-ইনডেক্সড) — র‍্যান্ডম পিভট, expected O(n)।"""
    a = list(a)
    lo, hi = 0, len(a) - 1
    while True:
        if lo == hi:
            return a[lo]
        pivot_idx = rng.randint(lo, hi)
        a[pivot_idx], a[hi] = a[hi], a[pivot_idx]
        pivot = a[hi]
        store = lo
        for i in range(lo, hi):
            if a[i] < pivot:
                a[store], a[i] = a[i], a[store]
                store += 1
        a[store], a[hi] = a[hi], a[store]
        if k == store:
            return a[store]
        elif k < store:
            hi = store - 1
        else:
            lo = store + 1

def find_median_quickselect(arr, rng):
    n = len(arr)
    if n % 2 == 1:
        return quickselect(arr, n // 2, rng)
    lower = quickselect(arr, n // 2 - 1, rng)
    upper = quickselect(arr, n // 2, rng)
    return (lower + upper) / 2

def true_median(arr):
    s = sorted(arr)
    n = len(s)
    if n % 2 == 1:
        return s[n // 2]
    return (s[n // 2 - 1] + s[n // 2]) / 2

all_correct = True
for seed in [1, 2, 3, 4, 5, 6, 7, 8]:
    rng = random.Random(seed)
    size = rng.randint(1, 60)
    arr = [rng.randint(0, 1000) for _ in range(size)]
    computed = find_median_quickselect(arr, rng)
    expected = true_median(arr)
    correct = (computed == expected)
    all_correct = all_correct and correct
    print(f"seed={seed} n={size}: quickselect মিডিয়ান={computed}, sorted() মিডিয়ান={expected}, মিলেছে? {correct}")

print(f"\nসব সিডে quickselect ও sorted()-ভিত্তিক মিডিয়ান মিলেছে: {all_correct}")
assert all_correct