PCA-র বাইরে — kernel PCA, ICA আর NMF

AI-এর ভিত্তি কোর্সের পর্ব ৬৫-এ PCA আপনি SVD দিয়ে চালিয়েছেন, explained variance আর পুনর্গঠন দেখেছেন; এই কোর্সের পর্ব ০৮-এ দেখেছেন PCA আসলে covariance matrix-এর eigenvector, আর পর্ব ০৯-এ truncated SVD দিয়ে low-rank approximation. এই পর্ব সেখান থেকে শুরু করে তিনটি প্রশ্নের উত্তর খোঁজে, যা PCA দিতে পারে না: গঠন বাঁকা হলে কী করবেন, কয়েকটি স্বাধীন উৎস মিশে গেলে কীভাবে আলাদা করবেন, আর ডেটা যদি শুধু অঋণাত্মক অংশ যোগ করে তৈরি হয়। চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python, scikit-learn 1.8.0); প্রথমবার scikit-learn নামতে প্রায় দশ সেকেন্ড লাগে।

১ · PCA কোথায় থামে — চারটি পদ্ধতি, চারটি প্রশ্ন

PCA-র উত্তর তিনটি শর্তে বাঁধা। এক, সে linear: নতুন অক্ষগুলো মূল feature-এর সরল ওজনযুক্ত যোগ, তাই বৃত্ত বা পেঁচানো গঠন খুলতে পারে না। দুই, সে শুধু দ্বিতীয় moment (variance আর covariance) দেখে, আর অক্ষগুলো পরস্পর লম্ব হতে বাধ্য; ফলে উপাদানগুলো uncorrelated হয়, কিন্তু স্বাধীন নয়। তিন, অক্ষে ঋণাত্মক ওজন থাকে: “চাল +০.৫, বই −০.৩” ধরনের উপাদান গণিতে ঠিক, কিন্তু কেনাকাটার ভাষায় পড়া কঠিন। নিচের ছবিতে চারটি পদ্ধতি একই ধরনের প্রশ্নের চারটি রূপ।

Probabilistic PCA আর factor analysis — এক অনুচ্ছেদে

PCA-কে একটি সম্ভাবনার মডেল হিসেবেও লেখা যায়। ধরুন প্রতিটি ডেটা বিন্দু x = Wz + μ + ε, যেখানে z হলো k মাত্রার গোপন কারণ (latent factor), z ~ N(০, I); W একটি d × k matrix; μ গড়; আর ε হলো noise, ε ~ N(০, σ²I). Michael Tipping আর Christopher Bishop ১৯৯৯ সালে দেখান, এই মডেলের maximum likelihood (পর্ব ১৮) সমাধানে W-এর কলামগুলো ঠিক প্রথম k টি principal component-এর জগৎ জুড়ে থাকে। এর নাম probabilistic PCA. লাভ: একটি likelihood পাওয়া যায় (model তুলনা আর অনুপস্থিত মান সামলানো সহজ), আর x-এর বণ্টন একটি পর্ব ১৯-এর multivariate Gaussian, যার covariance WWᵀ + σ²I. Noise-কে প্রতিটি feature-এ আলাদা হতে দিলে, অর্থাৎ σ²I-এর জায়গায় একটি diagonal matrix Ψ, পাওয়া যায় পুরোনো factor analysis. অংশ ৮-এ দুটিকে কাজে দেখবেন।

২ · Kernel PCA — Gram matrix কেন্দ্রে এনে eigenvector

পর্ব ২৯-এ দেখেছেন, যে algorithm শুধু dot product ব্যবহার করে, তাতে x·z-এর জায়গায় একটি kernel k(x, z) = φ(x)·φ(z) বসালে φ না বানিয়েই উঁচু মাত্রায় কাজ করা যায়। Bernhard Schölkopf, Alexander Smola আর Klaus-Robert Müller ১৯৯৮ সালে (Neural Computation) দেখান, PCA-ও এমন একটি algorithm. কারণটি সহজ: n টি বিন্দুর ডেটা matrix X (গড় শূন্য) হলে covariance XᵀX/n-এর eigenvector যে দিক দেয়, Gram matrix K = XXᵀ-এর eigenvector সেই একই দিকে প্রতিটি বিন্দুর অবস্থান দেয় (পর্ব ০৯-এর SVD-এর দুই পাশ)। আর K-এর প্রতিটি ঘর একটি dot product. তাই K_ij = k(x_i, x_j) বসালেই φ-এর জগতে PCA, যার নাম kernel PCA.

একটি বাধা আছে। PCA-র আগে ডেটার গড় শূন্য করতে হয়, কিন্তু φ(x) আমাদের হাতে নেই, তাই φ-এর গড় বাদ দেওয়া সরাসরি সম্ভব নয়। সমাধান হলো K-কেই কেন্দ্রে আনা (kernel centring). 1_n দিয়ে বোঝাই n × n matrix, যার প্রতিটি ঘর ১/n. তাহলে

K̃ = K − 1_nK − K1_n + 1_nK1_n

ঘর ধরে পড়লে: K̃_ij = K_ij − (i-তম সারির গড়) − (j-তম কলামের গড়) + (সব ঘরের গড়)। এর নাম double centring. ছোট একটি উদাহরণে দেখুন, যেখানে উত্তর আগে থেকেই জানা: linear kernel-এ K̃ হওয়া উচিত গড় বাদ দেওয়া ডেটার Gram matrix.

বাকি কাজ সাধারণ eigen-decomposition (পর্ব ০৮)। K̃α = λα থেকে eigenvalue বড় থেকে ছোট ক্রমে সাজান। প্রতিটি α-কে একক দৈর্ঘ্যের রাখলে train বিন্দুগুলোর projection হলো √λ · α. K̃ হলো পর্ব ১০-এর PSD matrix, তাই সব λ ≥ ০ (গণনার গোলমালে খুব ছোট ঋণাত্মক আসতে পারে)। Covariance-এর eigenvalue হলো λ/n. নতুন একটি বিন্দু x-এর জন্য লাগে k(x, x_i)-এর একটি সারি, যেটিও train-এর গড় দিয়ে একইভাবে কেন্দ্রে আনতে হয়; scikit-learn-এর KernelPCA.transform এটি নিজেই করে।

৩ · সমকেন্দ্রিক বৃত্তে kernel PCA, হাতে আর scikit-learn-এ

নিচের ব্লক scikit-learn-এর make_circles দিয়ে ২০০টি বিন্দুর দুটি সমকেন্দ্রিক বৃত্ত বানায় (seed দেওয়া)। তারপর হাতে RBF Gram matrix, double centring আর scipy.linalg.eigh দিয়ে kernel PCA করে, KernelPCA-এর সঙ্গে মেলায়, আর linear PCA-র সঙ্গে তুলনা করে। তুলনার মাপ: প্রথম উপাদানের অক্ষে একটিমাত্র কাটা দাগ দিয়ে দুই বৃত্ত সবচেয়ে ভালো কতটা আলাদা হয়। লক্ষ করুন, PCA নিজে লেবেল দেখে না; লেবেল শুধু ফল মাপতে।

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh
from sklearn.datasets import make_circles
from sklearn.decomposition import PCA, KernelPCA

X, y = make_circles(n_samples=200, factor=0.3, noise=0.05, random_state=41)
gamma, n = 2.0, len(X)
sq = ((X[:, None, :] - X[None, :, :]) ** 2).sum(axis=2)   # সব জোড়ার দূরত্বের বর্গ
K = np.exp(-gamma * sq)                                   # RBF Gram matrix, 200 x 200
J = np.full((n, n), 1 / n)
Kc = K - J @ K - K @ J + J @ K @ J                        # double centring
lam, A = eigh(Kc)                                         # ছোট থেকে বড় ক্রমে
lam, A = lam[::-1][:2], A[:, ::-1][:, :2]                 # সবচেয়ে বড় দুটি
Z_mine = A * np.sqrt(lam)                                 # projection = √λ · α
kp = KernelPCA(n_components=2, kernel="rbf", gamma=gamma).fit(X)
Z_skl = kp.transform(X)
sign = np.sign((Z_mine * Z_skl).sum(axis=0))              # eigenvector-এর চিহ্ন ইচ্ছামতো
print("my eigenvalues     :", np.round(lam, 3))
print("sklearn eigenvalues:", np.round(kp.eigenvalues_, 3))
print("sign flips:", sign, " max |diff| =", f"{np.abs(Z_mine * sign - Z_skl).max():.1e}")

def split_acc(z):                                         # এক অক্ষে একটি কাটা দাগে সেরা accuracy
    return max(max(((z > c) == y).mean(), ((z <= c) == y).mean()) for c in z)

Z_lin = PCA(n_components=2).fit_transform(X)
for name, Z in [("linear PCA", Z_lin), ("kernel PCA", Z_mine)]:
    print(f"{name}: one cut on component 1 -> accuracy {split_acc(Z[:, 0]):.3f}")

fig, ax = plt.subplots(1, 3, figsize=(11, 3.6))
for a, Z, t in zip(ax, [X, Z_lin, Z_mine], ["original", "linear PCA", "RBF kernel PCA"]):
    a.scatter(Z[:, 0], Z[:, 1], c=y, cmap="coolwarm", s=10)
    a.set(title=t, xlabel="component 1", ylabel="component 2")
fig.subplots_adjust(wspace=0.35, bottom=0.17)

my eigenvalues     : [29.661 24.46 ]
sklearn eigenvalues: [29.661 24.46 ]
sign flips: [-1.  1.]  max |diff| = 1.1e-14
linear PCA: one cut on component 1 -> accuracy 0.690
kernel PCA: one cut on component 1 -> accuracy 1.000

৪ · ICA — মেশানো গলা থেকে আলাদা বক্তা

Independent component analysis (ICA)-র মডেল খুব সরল: x = As. এখানে s হলো কয়েকটি অদেখা উৎসের মান (যেমন তিন বাবুর্চির গলা), A একটি অজানা মিশ্রণের matrix (কোন ফোন কার গলা কতটা শোনে), আর x হলো যা আমরা মাপি। লক্ষ্য: শুধু x দেখে A আর s দুটিই বের করা। A বা s কোনোটি না জেনে উৎস আলাদা করার এই কাজের নাম blind source separation. একমাত্র অনুমান: উৎসগুলো পরিসংখ্যানে স্বাধীন, আর (বড়জোর একটি ছাড়া) non-Gaussian.

প্রথম ধাপ whitening: x-কে এমনভাবে রূপান্তর করা যাতে covariance একক matrix হয় (পর্ব ১৯-এর হেলানো উপবৃত্ত গোল করা)। এটি PCA দিয়েই হয়, আর এতে কাজের অর্ধেক শেষ: whitening-এর পর বাকি থাকে শুধু একটি অজানা ঘূর্ণন। কিন্তু কোন ঘূর্ণন? দ্বিতীয় moment আর কিছু বলতে পারে না, কারণ whitened ডেটার যেকোনো ঘূর্ণনের covariance-ও একক matrix. এখানে ঢোকে উঁচু moment, বিশেষ করে kurtosis.

Kurtosis সহজ, কিন্তু কয়েকটি বড় মানের প্রতি খুব সংবেদনশীল। আরও শক্ত মাপ negentropy: একই variance-এর Gaussian-এর entropy থেকে আসল বণ্টনের entropy বাদ দেওয়া (AI-এর ভিত্তি কোর্সের পর্ব ৩৮-এ entropy). একই variance-এ Gaussian-এর entropy সবচেয়ে বেশি, তাই negentropy সবসময় ≥ ০, আর শূন্য শুধু Gaussian-এ। Aapo Hyvärinen আর Erkki Oja ১৯৯৭ সালে (Neural Computation) FastICA প্রকাশ করেন, যা negentropy-র একটি আনুমানিক রূপ (scikit-learn-এ ডিফল্ট fun="logcosh", অর্থাৎ G(u) = log cosh u) একটি fixed-point পুনরাবৃত্তিতে সর্বোচ্চ করে: w ← E[x g(wᵀx)] − E[g′(wᵀx)] w, তারপর w-কে একক দৈর্ঘ্যে আনা আর অন্য উপাদানের থেকে আলাদা (decorrelate) রাখা। এখানে g হলো G-এর অন্তরজ, অর্থাৎ tanh.

৫ · তিনটি সংকেত মিশিয়ে আবার আলাদা — PCA বনাম FastICA

নিচের ব্লকে তিনটি উৎস: একটি sine তরঙ্গ, একটি square তরঙ্গ আর একটি sawtooth (করাতের দাঁতের মতো) তরঙ্গ, সঙ্গে সামান্য noise. একটি 3 × 3 matrix A দিয়ে মিশিয়ে তিনটি “মাইক্রোফোন”-এর সংকেত বানানো হয়েছে। তারপর PCA (whitening-সহ) আর FastICA দুটিই চালিয়ে দেখা হয়েছে, পাওয়া উপাদান আর আসল উৎসের correlation কত। scipy.stats.kurtosis ডিফল্টে excess kurtosis দেয় (Gaussian-এ ০)।

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import kurtosis
from sklearn.decomposition import PCA, FastICA

rng = np.random.default_rng(41)
t = np.linspace(0, 8, 2000)
S = np.column_stack([np.sin(2 * t),                       # sine
                     np.sign(np.sin(3 * t)),              # square
                     2 * ((0.8 * t) % 1) - 1])            # sawtooth
S += 0.02 * rng.normal(size=S.shape)                      # সামান্য noise
A = np.array([[1.0, 1.0, 1.0], [0.5, 2.0, 1.0], [1.5, 1.0, 2.0]])
X = S @ A.T                                               # তিনটি “মাইক্রোফোন”
print("excess kurtosis of sources:", np.round(kurtosis(S), 2))
print("excess kurtosis of mixes  :", np.round(kurtosis(X), 2))

ica = FastICA(n_components=3, whiten="unit-variance", random_state=4)
S_ica = ica.fit_transform(X)
S_pca = PCA(n_components=3, whiten=True).fit_transform(X)
np.set_printoptions(suppress=True)
for name, E in [("PCA", S_pca), ("ICA", S_ica)]:
    C = np.corrcoef(S.T, E.T)[:3, 3:]                     # সারি = আসল উৎস, কলাম = পাওয়া উপাদান
    print(f"\n{name}: corr(true source, component)\n{np.round(C, 2)}")
    print("best |corr| per source:", np.round(np.abs(C).max(axis=1), 3))
print("\nstd: true", np.round(S.std(axis=0), 2), "| ICA", np.round(S_ica.std(axis=0), 2))

fig, ax = plt.subplots(3, 3, figsize=(10, 5.5), sharex=True)
for j in range(3):
    ax[j, 0].plot(t[:600], S[:600, j]); ax[j, 1].plot(t[:600], X[:600, j])
    ax[j, 2].plot(t[:600], S_ica[:600, j], color="tab:green")
for a, ttl in zip(ax[0], ["true sources", "mixtures (observed)", "FastICA output"]):
    a.set_title(ttl)
fig.subplots_adjust(hspace=0.3, wspace=0.25)

excess kurtosis of sources: [-1.36 -1.99 -1.21]
excess kurtosis of mixes  : [-1.17 -1.68 -1.21]

PCA: corr(true source, component)
[[ 0.46 -0.57  0.67]
 [ 0.8   0.59  0.08]
 [ 0.43 -0.55 -0.72]]
best |corr| per source: [0.674 0.804 0.716]

ICA: corr(true source, component)
[[ 0.06 -1.   -0.05]
 [ 0.03 -0.03 -1.  ]
 [ 1.    0.03  0.06]]
best |corr| per source: [0.997 0.999 0.998]

std: true [0.69 1.   0.58] | ICA [1. 1. 1.]

scikit-learn 1.8.0-এ whiten দুটি string নেয়: "unit-variance" (ডিফল্ট) বা "arbitrary-variance", অথবা False. পুরোনো টিউটোরিয়ালের whiten=True এখন InvalidParameterError দেয় (ডেস্কটপে 1.8.0-এ যাচাই করা)। বাস্তবে ICA-র বড় সীমা হলো মডেলটি নিজে: ঘরের সত্যিকারের রেকর্ডিংয়ে শব্দ দেয়ালে প্রতিফলিত হয়ে দেরিতে পৌঁছায়, তাই মিশ্রণ আর একটি সরল matrix নয় (convolutive mixing); আর মাইক্রোফোনের চেয়ে উৎস বেশি হলে x = As উল্টানো যায় না। EEG আর fMRI-এর মতো ক্ষেত্রে, যেখানে মিশ্রণ প্রায় তাৎক্ষণিক আর রৈখিক, ICA সবচেয়ে ভালো কাজ করে।

৬ · NMF — শুধু যোগ করে অংশ থেকে সম্পূর্ণ

Non-negative matrix factorization (NMF) একটি অঋণাত্মক n × d matrix V-কে দুটি অঋণাত্মক matrix-এর গুণফল হিসেবে আন্দাজ করে: V ≈ WH, যেখানে W হলো n × k (প্রতিটি ক্রেতার প্রতিটি অংশে কতটা ওজন) আর H হলো k × d (প্রতিটি অংশে কোন বিভাগ কতটা)। এটি পর্ব ০৯-এর matrix factorization, শুধু একটি বাড়তি শর্তে: W ≥ ০ আর H ≥ ০, ঘরে ঘরে। লক্ষ্য Frobenius norm-এর বর্গ কমানো:

সর্বনিম্ন করুন ‖V − WH‖_F² = ∑_i,j (V_ij − (WH)_ij)²,  শর্ত W ≥ ০, H ≥ ০

যেহেতু কোনো বিয়োগ নেই, একটি অংশ আরেকটিকে কাটাকাটি করতে পারে না; ফলে অংশগুলো সাধারণত sparse আর পড়ার মতো হয়। একে parts-based representation বলে। Pentti Paatero আর Unto Tapper ১৯৯৪ সালে “positive matrix factorization” নামে বায়ুদূষণের উৎস খোঁজায় এই ধারণা ব্যবহার করেন; Daniel Lee আর Sebastian Seung ১৯৯৯ সালে Nature-এ মুখের ছবি থেকে চোখ, নাক, ঠোঁটের মতো অংশ শেখার উদাহরণ দিয়ে একে বিখ্যাত করেন, আর ২০০১ সালে সহজ multiplicative update নিয়ম প্রকাশ করেন:

H ← H ∘ (WᵀV) ⊘ (WᵀWH),   W ← W ∘ (VHᵀ) ⊘ (WHHᵀ)

এখানে ∘ মানে ঘরে ঘরে গুণ, ⊘ মানে ঘরে ঘরে ভাগ। নিয়মটি আসলে gradient descent, যার প্রতিটি ঘরের ধাপের মাপ এমনভাবে বাছা যে বিয়োগ চিহ্ন গুণ-ভাগে বদলে যায়। শুরুর মান ধনাত্মক হলে মান কখনো ঋণাত্মক হয় না, আর Lee ও Seung প্রমাণ করেন প্রতিটি ধাপে ত্রুটি বাড়ে না। নিচে একটি 2 × 2 উদাহরণে একটি পুরো ধাপ:

৭ · বাজারঘরের কেনাকাটা ভেঙে প্রোফাইল

এই অংশের ডেটা কাল্পনিক, seed দিয়ে বানানো। বাজারঘরের ৮০ জন ক্রেতা, আটটি পণ্য-বিভাগ (চাল-ডাল, মাছ-মাংস, সবজি, মসলা, শিশুর পণ্য, প্রসাধনী, গ্যাজেট, বই), আর প্রতিটি ঘরে তিন মাসে সেই বিভাগে কতবার কেনা হয়েছে। ভেতরে তিনটি লুকানো অভ্যাস আছে, যা আমরা জানি কিন্তু algorithm জানে না। ব্লক ৩ আর ৪ এই module থেকে ডেটা নেয়।

# bazar_basket.py — বাজারঘরের কাল্পনিক ক্রেতা × পণ্য-বিভাগ কেনাকাটা (seed দিয়ে বানানো)
import numpy as np

CATS = ["rice_dal", "fish_meat", "vegetable", "spice", "baby", "cosmetic", "gadget", "book"]

def make_basket(n=80, seed=41):
    rng = np.random.default_rng(seed)
    H = np.array([[6, 4, 5, 3, 0, 0.5, 0, 0],      # সংসারের বাজার
                  [0.5, 0, 0.5, 0, 6, 3, 0, 0.5],  # নতুন বাবা-মা
                  [0, 0, 0, 0, 0, 1, 4, 5]])       # ছাত্র-ছাত্রী, গ্যাজেট
    W = rng.gamma(0.6, 1.0, size=(n, 3))           # প্রতিজন ক্রেতায় তিন অভ্যাসের মিশ্রণ
    V = rng.poisson(W @ H).astype(float)           # তিন মাসে প্রতিটি বিভাগে কেনার সংখ্যা
    return V, CATS

নিচের ব্লক multiplicative update হাতে লেখে (১,০০০ ধাপ, কয়েকটি ধাপে ত্রুটি ছাপে), তারপর scikit-learn-এর NMF-এর সঙ্গে ত্রুটি মেলায়। scikit-learn-এর reconstruction_err_ হলো ‖V − WH‖_F, বর্গ নয়। init="nndsvda" একটি SVD-ভিত্তিক শুরু, আর ডিফল্ট solver "cd" (coordinate descent); max_iter যথেষ্ট রাখায় convergence-এর কোনো সতর্কবার্তা আসে না।

import numpy as np
from sklearn.decomposition import NMF
from bazar_basket import make_basket

V, cats = make_basket()                                   # ৮০ ক্রেতা × ৮ বিভাগ
print("V shape:", V.shape, " zeros:", f"{(V == 0).mean():.0%}")

def nmf_mu(V, k, iters=1000, seed=0, eps=1e-10):          # Lee–Seung multiplicative update
    rng = np.random.default_rng(seed)
    W = rng.random((V.shape[0], k)); H = rng.random((k, V.shape[1]))
    for it in range(iters):
        H *= (W.T @ V) / (W.T @ W @ H + eps)
        W *= (V @ H.T) / (W @ H @ H.T + eps)
        if it in (0, 9, 99, 999):
            print(f"iter {it + 1:4d}  ||V - WH||_F = {np.linalg.norm(V - W @ H):.4f}")
    return W, H

W, H = nmf_mu(V, k=3)
skl = NMF(n_components=3, init="nndsvda", solver="cd", max_iter=1000, random_state=0)
W2 = skl.fit_transform(V)
print(f"sklearn NMF (cd)          error = {skl.reconstruction_err_:.4f}, iters {skl.n_iter_}")
print("all entries of W and H >= 0:", bool((W >= 0).all() and (H >= 0).all()))

Hn = H / H.sum(axis=1, keepdims=True)                     # প্রতিটি অংশ একটি কেনাকাটার প্রোফাইল
for r, row in enumerate(Hn):
    top = np.argsort(row)[::-1][:3]
    print(f"part {r}: " + ", ".join(f"{cats[j]} {row[j]:.2f}" for j in top))
c = 7
print(f"\ncustomer {c}: buys", V[c].astype(int), "\n  weights W =", W[c].round(2))

V shape: (80, 8)  zeros: 30%
iter    1  ||V - WH||_F = 72.8120
iter   10  ||V - WH||_F = 34.6105
iter  100  ||V - WH||_F = 32.4645
iter 1000  ||V - WH||_F = 32.4475
sklearn NMF (cd)          error = 32.4476, iters 33
all entries of W and H >= 0: True
part 0: book 0.49, gadget 0.36, cosmetic 0.12
part 1: baby 0.59, cosmetic 0.34, vegetable 0.04
part 2: rice_dal 0.33, vegetable 0.28, fish_meat 0.21

customer 7: buys [ 8  5  6  4  0  0 12 15]
  weights W = [2.33 0.   1.38]

৮ · PCA, NMF আর factor analysis পাশাপাশি — আর কখন কোনটি

একই কেনাকাটার ডেটায় তিনটি পদ্ধতি চালিয়ে উপাদানগুলো পাশাপাশি দেখা যাক, আর k বাড়ালে পুনর্গঠনের ত্রুটি (reconstruction error) কীভাবে নামে। তিন রকম ত্রুটি তুলনা করা হয়েছে: SVD-এর সেরা rank-k approximation (গড় বাদ না দিয়ে), PCA (গড় বাদ দিয়ে k টি উপাদান, তারপর গড় ফেরত), আর NMF. Factor analysis চালানো হয়েছে varimax rotation দিয়ে, যা factor-গুলোকে এমনভাবে ঘোরায় যাতে প্রতিটি factor অল্প কয়েকটি feature-এ বড় ওজন পায়।

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import NMF, PCA, FactorAnalysis
from bazar_basket import make_basket

V, cats = make_basket()
err = lambda R: np.linalg.norm(V - R)                     # Frobenius ত্রুটি
s = np.linalg.svd(V, compute_uv=False)                    # singular value, বড় থেকে ছোট
print(" k  SVD(best rank k)  PCA(mean + k)  NMF")
for k in range(1, 7):
    p = PCA(n_components=k).fit(V)
    m = NMF(n_components=k, init="nndsvda", max_iter=2000, random_state=0)
    Wk = m.fit_transform(V)
    print(f"{k:2d}  {np.sqrt((s[k:] ** 2).sum()):14.2f}  {err(p.inverse_transform(p.transform(V))):13.2f}"
          f"  {err(Wk @ m.components_):6.2f}")

pca = PCA(n_components=3).fit(V)
nmf = NMF(n_components=3, init="nndsvda", max_iter=2000, random_state=0).fit(V)
fa = FactorAnalysis(n_components=3, rotation="varimax", random_state=0).fit(V)
print("\nnegative entries: PCA", (pca.components_ < 0).sum(), "of 24 | NMF",
      (nmf.components_ < 0).sum(), "| FA", (fa.components_ < 0).sum())
print("PCA explained variance ratio:", pca.explained_variance_ratio_.round(3))
print("FA noise variance per category:", fa.noise_variance_.round(2))

fig, ax = plt.subplots(1, 3, figsize=(11, 3.4))
for a, M, t in zip(ax, [pca.components_, nmf.components_, fa.components_],
                   ["PCA components", "NMF parts", "Factor analysis (varimax)"]):
    M = M / np.abs(M).max(axis=1, keepdims=True)          # প্রতিটি সারি নিজের সর্বোচ্চে ভাগ
    a.imshow(M, cmap="RdBu_r", vmin=-1, vmax=1, aspect="auto")
    a.set(title=t, yticks=range(3), xticks=range(8))
    a.set_xticklabels(cats, rotation=60, fontsize=8)
fig.subplots_adjust(wspace=0.3, bottom=0.32)

 k  SVD(best rank k)  PCA(mean + k)  NMF
 1           73.07          68.35   73.07
 2           52.80          52.70   53.02
 3           32.14          32.08   32.45
 4           26.60          26.59   27.00
 5           21.15          21.13   22.92
 6           15.97          15.96   17.28

negative entries: PCA 8 of 24 | NMF 0 | FA 9
PCA explained variance ratio: [0.428 0.232 0.214]
FA noise variance per category: [4.72 3.07 2.13 0.92 1.67 3.66 2.28 1.75]

কোন প্রশ্নে কোন পদ্ধতি

শুধু ছবি আঁকার জন্য, যেখানে কাছের প্রতিবেশী ঠিক রাখাই লক্ষ্য, পর্ব ৪২ (t-SNE আর UMAP — manifold learning আর তার ফাঁদ) দেখবেন। আর মনে রাখবেন, NMF দিয়ে লেখার বিষয় খোঁজা পর্ব ৪৪ (Topic modeling — LDA আর NMF দিয়ে বাংলা রিভিউয়ের বিষয় খোঁজা)-এর বিষয়।

৯ · AI-এর যুগে — dimensionality reduction-এর কোড যাচাই

AI assistant-কে “এই ডেটায় ICA বা NMF চালাও” বললে সে সাধারণত চলনসই কোড দেয়, কিন্তু ভুল থাকে পুরোনো API, উপাদানের অর্থ আর মাপার পদ্ধতিতে। পরিচিত ভুলগুলো:

AI-এর দেওয়া dimensionality reduction কোড হাতে পেলে — চেকলিস্ট

• প্রতিটি argument নিজের scikit-learn version-এর docs-এ মিলিয়েছি, আর কোড চালিয়ে কোনো warning চাপা পড়েনি
• Shape আর কয়েকটি সারি দেখে নিশ্চিত হয়েছি কোনটি নমুনা আর কোনটি feature
• NMF-এর আগে ডেটা অঋণাত্মক; ICA-র আগে উৎসগুলো non-Gaussian হওয়ার কারণ আছে
• Model-এর আগে কমানো হলে fit শুধু train-এ, Pipeline-এর ভেতরে, আর test তালাবন্ধ
• কয়েকটি seed আর k-এ উপাদানের স্থিরতা দেখেছি; একটি run-এর অংশকে সত্য ধরিনি
• সূত্র (যেমন double centring বা multiplicative update) একটি ছোট উদাহরণে হাতে মিলিয়েছি
• আসল গ্রাহকের ব্যক্তিগত ডেটা AI-তে পেস্ট করিনি; নমুনা বা কাল্পনিক ডেটা ব্যবহার করেছি