Clustering আর PCA — label ছাড়া গঠন খোঁজা

এ পর্যন্ত প্রায় সব model-এর হাতে উত্তর (label) ছিল; এই পর্বে উত্তর নেই, শুধু ডেটা। চালানো যায় এমন ব্লকগুলো পাতাতেই চলে (ব্রাউজারের Python: NumPy 2.4.6, SciPy 1.18.0, scikit-learn 1.8.0, matplotlib 3.10.8). scikit-learn প্রথমবার নামতে প্রায় ১০ সেকেন্ড লাগে; পরের run দ্রুত।

১ · label ছাড়া শেখা — বাজারঘরের গ্রাহকদের দল

পর্ব ৫৫ (যন্ত্র কীভাবে শেখে)-এ দেখেছেন, unsupervised learning-এ শুধু X থাকে, y থাকে না; model ডেটার ভেতরের গঠন খোঁজে। দুটি প্রধান কাজ: clustering (একই রকম সারিকে এক দলে রাখা) আর dimensionality reduction (অনেক কলামকে কম কলামে আনা, যতটা সম্ভব তথ্য রেখে)। এই পর্বে দুটিই, হাতে-কলমে।

বাজারঘর (এই কোর্সের কাল্পনিক অনলাইন দোকান) তার গ্রাহকদের জন্য আলাদা বার্তা পাঠাতে চায়: নিয়মিত ক্রেতাকে ধন্যবাদ, ঘুমিয়ে পড়া ক্রেতাকে ফিরে আসার অনুরোধ। প্রতিটি গ্রাহক তিনটি সংখ্যায়, direct marketing-এ বহু দশক ধরে প্রচলিত RFM ধারণা:

• Recency (R): শেষ কেনার পর কত দিন গেছে। কম মানে সম্প্রতি কিনেছেন।
• Frequency (F): বছরে কয়টি অর্ডার।
• Monetary (M): বছরে মোট কত টাকার কেনাকাটা।

আমাদের ডেটা কাল্পনিক, seed দিয়ে বানানো: ৩০০ জন গ্রাহক, চারটি লুকানো দল থেকে। (১) নিয়মিত ছোট কেনাকাটা: ৯০ জন, R প্রায় ১২ দিন, F প্রায় ১৪, M প্রায় ৯০০০ টাকা; (২) কম কিন্তু বড় কেনাকাটা: ৫০ জন, R ৪০, F ৪, M ২৪০০০; (৩) ঘুমিয়ে পড়া: ১০০ জন, R ২১০, F ২, M ২৫০০; (৪) মাঝে মাঝে ছোট কেনাকাটা: ৬০ জন, R ২৫, F ৩, M ৩৫০০. দলগুলো আমরা জানি, algorithm জানে না; শেষে মিলিয়ে দেখব। বাস্তবের ডেটায় এমন “লুকানো সঠিক দল” প্রায়ই থাকে না।

২ · k-means — দুই ধাপের loop, হাতে এক iteration

k-means-এর লক্ষ্য একটি সংখ্যা ছোট করা। প্রতিটি দলের একটি কেন্দ্র (centroid) থাকে, আর প্রতিটি বিন্দু তার নিজের দলের কেন্দ্র থেকে যত দূরে, সেই Euclidean distance-এর বর্গের যোগফলই inertia (within-cluster sum of squares):

J = ∑_i ‖x_i − c_দল(i)‖²

এখানে x_i হলো i-তম বিন্দু, c_দল(i) তার দলের কেন্দ্র, ‖·‖² বর্গ দূরত্ব। সব সম্ভাব্য ভাগ চেষ্টা করা অসম্ভব রকম ব্যয়বহুল, তাই k-means দুটি সহজ ধাপ পালাক্রমে চালায় (Lloyd-এর algorithm):

1. Assign: কেন্দ্র স্থির রেখে প্রতিটি বিন্দুকে সবচেয়ে কাছের কেন্দ্রের দলে দেওয়া। এতে J বাড়ে না।
2. Update: দল স্থির রেখে প্রতিটি কেন্দ্রকে তার দলের বিন্দুগুলোর গড়ে সরানো। এতেও J বাড়ে না, কারণ বর্গ দূরত্বের যোগফল গড়েই সবচেয়ে ছোট (পর্ব ৩৭)।

দুই ধাপেই J কমে বা সমান থাকে, আর সম্ভাব্য ভাগের সংখ্যা সসীম, তাই loop একসময় থামে: কোনো বিন্দুর দল আর বদলায় না। কিন্তু থামার জায়গাটি সবচেয়ে ছোট J না-ও হতে পারে; এটি একটি local minimum, ঠিক পর্ব ৪৪-এর hill climbing-এর মতো।

ছয়টি বিন্দু আর k = ২ দিয়ে হাতে, যেমন বিশ্ববিদ্যালয়, চাকরি বা ইন্টারভিউর পরীক্ষায় আসে। শুরুর কেন্দ্র ইচ্ছে করে খারাপ বাছা: c1 = P1 = (১, ১) আর c2 = P3 = (১, ৩), দুটিই বাঁ দিকের দলে।

এখানে খারাপ শুরু থেকেও ঠিক ভাগ এসেছে, কারণ দল দুটি অনেক দূরে; অংশ ৩-এ দেখব সব সময় এমন হয় না।

৩ · NumPy-তে k-means — k-means++, restart আর scaling

কেন scaling আগে

k-means দূরত্ব দিয়ে চলে, তাই পর্ব ৩৩ আর পর্ব ৫৭-এর শিক্ষা এখানেও খাটে: কাঁচা সংখ্যায় monetary (হাজার হাজার টাকা) recency (শত দিন) আর frequency (কয়েকটি অর্ডার)-কে ডুবিয়ে দেয়। তখন “দল” মানে কার্যত শুধু খরচের ধাপ। তাই আগে standardization: প্রতিটি কলাম থেকে গড় বাদ দিয়ে standard deviation দিয়ে ভাগ। টাকার মতো ডান-দিকে-লম্বা কলামে log রূপান্তরও প্রচলিত; scaling-এর নানা পথ পর্ব ৬৩ (Feature engineering)-এর বিষয়।

k-means++: কেন্দ্র ছড়িয়ে শুরু

অংশ ২-এর মতো দুটি কেন্দ্র এক দলে পড়লে কষ্ট বাড়ে। k-means++ (David Arthur আর Sergei Vassilvitskii, ২০০৭) শুরুর কেন্দ্র এভাবে বাছে: প্রথমটি এলোমেলো একটি বিন্দু; তারপর প্রতিটি বিন্দুর সবচেয়ে কাছের বাছা কেন্দ্র থেকে বর্গ দূরত্ব D(x)² মাপা হয়, আর পরের কেন্দ্র বাছা হয় D(x)²-এর অনুপাতে সম্ভাবনা দিয়ে। দূরের বিন্দুর কেন্দ্র হওয়ার সম্ভাবনা বেশি, তাই কেন্দ্রগুলো ছড়িয়ে পড়ে। তাঁরা দেখান, এই শুরুতেই প্রত্যাশিত inertia সেরা সম্ভবের O(log k) গুণের মধ্যে থাকে; এটি গড়ের নিশ্চয়তা, একটি run তবু খারাপ হতে পারে।

Restart: local minimum-এর ওষুধ

পর্ব ৪৪-এর random-restart hill climbing-এর মতো: কয়েকবার (এখানে ১০) আলাদা শুরু থেকে চালিয়ে সবচেয়ে কম inertia-র run রাখা। scikit-learn-এ এর নাম n_init.

নিচের ব্লক ডেটা বানায়, standardize করে, k-means++ আর সরল এলোমেলো শুরু দুটোই লেখে, তারপর (১) k = ৪-এ দুই রকম শুরুর ১০টি করে run-এর inertia, (২) k = ২ থেকে ৭-এ ১০ restart-এর সেরা inertia আর silhouette score, (৩) k = ৪-এর দলগুলোর আসল এককে গড় ছাপে। scikit-learn-এর তৈরি KMeans এখানে ব্যবহার করা হয়নি, কারণ ব্রাউজারের Python-এ এটি threadpoolctl থেকে একটি অপ্রাসঙ্গিক deprecation warning ছাপে; নিজে লেখা k-means-এ ভেতরটাও দেখা যায়। এই ব্লক আর পরের দুটি একই namespace ভাগ করে (data-shared); পরেরটি আগে চালালে এটি আপনা থেকে আগে চলে।

import numpy as np
from sklearn.metrics import silhouette_score, adjusted_rand_score
rng = np.random.default_rng(42)
# কাল্পনিক গ্রাহক: recency (শেষ কেনার পর দিন), frequency (বছরে অর্ডার), monetary (বছরে টাকা)
mu = [[12, 14, 9000], [40, 4, 24000], [210, 2, 2500], [25, 3, 3500]]   # ৪টি লুকানো দল
sd = [[6, 3, 2000], [15, 1.5, 5000], [60, 1, 800], [10, 1, 900]]
X = np.vstack([rng.normal(m, s, (n, 3)) for m, s, n in zip(mu, sd, [90, 50, 100, 60])]).clip(1)
Z = (X - X.mean(0)) / X.std(0)                     # standardize: প্রতিটি কলাম গড় 0, σ 1
truth = np.repeat(np.arange(4), [90, 50, 100, 60])  # লুকানো দল; algorithm এটি দেখে না
def kmeans(Z, k, rng, plusplus=True, iters=100):
    if plusplus:                                   # k-means++: দূরের বিন্দু কেন্দ্র হওয়ার সম্ভাবনা বেশি
        C = Z[[rng.integers(len(Z))]]
        while len(C) < k:
            d2 = ((Z[:, None] - C[None]) ** 2).sum(-1).min(1)
            C = np.vstack([C, Z[rng.choice(len(Z), p=d2 / d2.sum())]])
    else:
        C = Z[rng.choice(len(Z), k, replace=False)]  # সরল: k-টি এলোমেলো বিন্দু
    for _ in range(iters):
        lab = ((Z[:, None] - C[None]) ** 2).sum(-1).argmin(1)          # ধাপ ১: assign
        newC = np.array([Z[lab == j].mean(0) if (lab == j).any() else C[j] for j in range(k)])
        if np.allclose(newC, C): break                                 # ধাপ ২: update
        C = newC
    return lab, C, ((Z - C[lab]) ** 2).sum()       # inertia: নিজের কেন্দ্র থেকে বর্গ দূরত্বের যোগ
rng = np.random.default_rng(0)
print("random init, ১০ বার:", np.sort([kmeans(Z, 4, rng, False)[2] for _ in range(10)]).round(1))
print("k-means++,  ১০ বার:", np.sort([kmeans(Z, 4, rng)[2] for _ in range(10)]).round(1))
print(" k  inertia  silhouette")
for k in range(2, 8):
    lab, C, inertia = min((kmeans(Z, k, rng) for _ in range(10)), key=lambda r: r[2])
    print(f"{k:2d} {inertia:8.1f}  {silhouette_score(Z, lab):.3f}")
lab4, C4, _ = min((kmeans(Z, 4, rng) for _ in range(10)), key=lambda r: r[2])
print("লুকানো দলের সঙ্গে মিল (ARI):", round(adjusted_rand_score(truth, lab4), 3))
for j in np.argsort(C4[:, 0]):                     # দল, recency-র ক্রমে, আসল এককে গড়
    r, f, m = X[lab4 == j].mean(0)
    print(f"দল {j}: {np.sum(lab4 == j):3d} জন | recency {r:5.0f} দিন | {f:4.1f} অর্ডার | {m:6.0f} টাকা")

random init, ১০ বার: [ 88.8  88.8  88.8  88.8  88.8  88.8  88.8 272.4 272.5 302.1]
k-means++,  ১০ বার: [ 88.8  88.8  88.8  88.8  88.8  88.8 226.1 226.1 226.1 230.8]
 k  inertia  silhouette
 2    473.9  0.483
 3    261.7  0.560
 4     88.8  0.667
 5     70.8  0.606
 6     54.6  0.543
 7     42.6  0.523
লুকানো দলের সঙ্গে মিল (ARI): 0.851
দল 3:  89 জন | recency    12 দিন | 14.0 অর্ডার |   8927 টাকা
দল 0:  78 জন | recency    40 দিন |  2.9 অর্ডার |   3490 টাকা
দল 1:  48 জন | recency    41 দিন |  3.9 অর্ডার |  24873 টাকা
দল 2:  85 জন | recency   222 দিন |  2.1 অর্ডার |   2440 টাকা

নিজের লেখা কোড ঠিক কি না, তা ডেস্কটপে scikit-learn-এর আসল KMeans-এর সঙ্গে মিলিয়ে দেখা হয়েছে: একই ডেটা আর একই Z, KMeans(n_clusters=4, n_init=10, random_state=0).

~/bazarghar-ai$ python kmeans_check.py
scikit-learn 1.9.1 | numpy 2.5.4
নিজের k-means inertia : 88.772
KMeans.inertia_        : 88.772
দুই ভাগের মিল (ARI)    : 1.0
KMeans কয়বার ঘুরল (n_iter_): 4

সত্যিকারের run (ডেস্কটপ, লেখার সময়ে Python 3.14.4, NumPy 2.5.4, scikit-learn 1.9.1). Inertia তিন দশমিক পর্যন্ত এক, আর দুই ভাগের ARI ১.০, অর্থাৎ প্রতিটি গ্রাহক একই দলে (নম্বর আলাদা হতে পারে)। সেরা run-এ মাত্র ৪টি iteration লেগেছে। একটি সূক্ষ্ম কথা: scikit-learn 1.4 থেকে n_init-এর default "auto", যা init="k-means++"-এর সঙ্গে মাত্র ১টি run চালায় (documentation অনুযায়ী)। উপরের ফলে আমাদের সরল k-means++-এর ১০টি run-এর ৪টি খারাপ ভাগে আটকেছে; তাই গুরুত্বপূর্ণ কাজে n_init নিজে বলে দেওয়াই নিরাপদ। (scikit-learn-এর k-means++ “greedy”: প্রতিটি কেন্দ্রের জন্য 2 + ⌊ln k⌋টি প্রার্থী থেকে সেরাটি নেয়, তাই তার একটি run আমাদের সরল রূপের চেয়ে কিছুটা নির্ভরযোগ্য।)

৪ · k কত? Elbow, silhouette আর দলের দায়িত্বশীল নাম

Inertia আর elbow

k বাড়ালে inertia প্রায় সবসময় কমে; k = n হলে প্রতিটি বিন্দু নিজেই কেন্দ্র, inertia ০. তাই “সবচেয়ে কম inertia” দিয়ে k বাছা যায় না। Elbow method দেখে কোথায় কমার গতি হঠাৎ ধীর হয়, হাতের কনুইয়ের মতো বাঁক। উপরের টেবিলে k = ৩ থেকে ৪-এ inertia ২৬১.৭ থেকে ৮৮.৮, প্রায় ১৭৩ কম; ৪ থেকে ৫-এ মাত্র ১৮. বাঁক k = ৪-এ। বাস্তবে বাঁক প্রায়ই এত স্পষ্ট নয়।

Silhouette: প্রতিটি বিন্দু কতটা নিজের দলে

Peter Rousseeuw ১৯৮৭ সালে প্রতিটি বিন্দু i-এর জন্য একটি মাপ প্রস্তাব করেন। a(i) = নিজের দলের অন্য বিন্দুগুলোর সঙ্গে গড় দূরত্ব; b(i) = সবচেয়ে কাছের অন্য দলের বিন্দুগুলোর সঙ্গে গড় দূরত্ব। তারপর

s(i) = (b(i) − a(i)) / max(a(i), b(i))

s(i)-এর মান −১ থেকে ১. ১-এর কাছে মানে নিজের দলের গায়ে লেগে আছে আর অন্য দল থেকে দূরে; ০-এর কাছে মানে দুই দলের সীমানায়; ঋণাত্মক মানে সম্ভবত ভুল দলে। সব বিন্দুর গড়ই silhouette score (sklearn.metrics.silhouette_score). আমাদের টেবিলে সর্বোচ্চ ০.৬৬৭, k = ৪-এ; elbow-এর সঙ্গে একমত। সতর্কতা: silhouette-ও গোল, ঘন দল পছন্দ করে, আর সব জোড়ার দূরত্ব লাগে বলে বড় ডেটায় ব্যয়বহুল।

দলের নাম দেওয়া — দায়িত্ব নিয়ে

ব্লকের শেষ চার লাইন দলগুলোর আসল এককে গড়। একজন বিশ্লেষক নাম দিতে পারেন: দল ৩ “নিয়মিত ক্রেতা” (R ১২ দিন, ১৪টি অর্ডার), দল ১ “বড় ঝুড়ির ক্রেতা” (প্রায় ২৫ হাজার টাকা), দল ০ “মাঝে মাঝে ছোট কেনাকাটা”, দল ২ “ঘুমিয়ে পড়া” (R ২২২ দিন)। কিন্তু মনে রাখুন:

• নাম একটি ব্যাখ্যা, তথ্য নয়। “বড় ঝুড়ির ক্রেতা” হয়তো অফিসের হয়ে কেনা মানুষ। আগে কয়েকটি সারি নিজে পড়ুন, শুধু গড় নয়।
• স্থিরতা যাচাই। অন্য seed, অন্য মাস বা ৮০% নমুনায় একই দল ফেরে কি না (দুই run-এর ARI). যে দল প্রতিবার বদলায়, তার নামে নীতি নয়।
• মানুষকে ছোট করে না এমন নাম। “অলস” বা “কম মূল্যের গ্রাহক” নয়; বর্ণনামূলক নাম। আর লিঙ্গ, ধর্ম বা এলাকার মতো সংবেদনশীল তথ্য দলের ভেতরে লুকিয়ে আছে কি না, তা আলাদা করে দেখুন (ভাবনার প্রশ্ন ৪)।
• দল ব্যবহারের ফল মাপুন। আলাদা বার্তায় কি সত্যিই ফেরত আসা বাড়ে? তা মাপে A/B test, silhouette নয়।

৫ · k-means যা ধরে নেয় — আর hierarchical clustering

k-means প্রতিটি বিন্দুকে সবচেয়ে কাছের কেন্দ্রে দেয়, তাই দুই কেন্দ্রের মাঝের সীমারেখা সবসময় একটি সরল রেখা (বা উচ্চ মাত্রায় সমতল), আর প্রতিটি দল একটি উত্তল (convex) এলাকা। এর ফলে কিছু অনুমান চুপচাপ ঢুকে পড়ে:

• দল মোটামুটি গোল: সব দিকে সমান ছড়ানো। লম্বাটে, তির্যক বা বাঁকা দল ভুলভাবে কাটা পড়ে।
• দলগুলোর আকার আর ছড়ানো কাছাকাছি: বড় ছড়ানো দলের পাশে ছোট ঘন দল থাকলে সীমারেখা ভুল জায়গায় পড়ে।
• প্রতিটি বিন্দু কোনো দলের: ছিটকে পড়া বিন্দুও দলে ঢোকে আর কেন্দ্রকে টেনে সরায়।
• k আগে থেকে জানা, আর দূরত্ব Euclidean.

অংশ ৬-এ দুটি বাঁকা “চাঁদ” (make_moons)-এ এই ব্যর্থতা দেখবেন। আগে আরেক ধারার পদ্ধতি।

Hierarchical clustering: জোড়া লাগাতে লাগাতে গাছ

Agglomerative hierarchical clustering শুরুতে প্রতিটি বিন্দুকে আলাদা দল ধরে, তারপর বারবার সবচেয়ে কাছের দুটি দল জোড়া লাগায়, যতক্ষণ না একটি দল থাকে (পর্ব ৫৫-এর ব্লকে একটি ছোট run দেখেছেন)। পুরো ইতিহাস আঁকা যায় একটি গাছে, যার নাম dendrogram: প্রতিটি জোড়া লাগার উচ্চতা হলো তখন দুই দলের দূরত্ব। k আগে বলতে হয় না; গাছটি কোনো উচ্চতায় আড়াআড়ি কাটলেই দল।

“দুই দলের দূরত্ব” কীভাবে মাপা হবে, তা ঠিক করে linkage. Single linkage-এ দুই দলের সবচেয়ে কাছের দুটি বিন্দুর দূরত্ব। নিচে এক মাত্রার পাঁচটি বিন্দু (যেমন পাঁচজন গ্রাহকের মাসিক অর্ডার) দিয়ে single linkage-এর গাছ বানিয়ে কাটা হলো।

উচ্চতাগুলো scipy 1.18-এর scipy.cluster.hierarchy.linkage দিয়ে ডেস্কটপে মাপা (ward-এর উচ্চতা scipy-র নিজস্ব সূত্রে, সাধারণ দূরত্ব নয়)। Ward-এর পদ্ধতি Joe H. Ward Jr.-এর ১৯৬৩ সালের প্রবন্ধ থেকে। Single linkage minimum spanning tree-র আত্মীয়: MST-র সবচেয়ে লম্বা k − ১টি edge কাটলে single-linkage-এর k দল মেলে (Gower ও Ross, ১৯৬৯)। সব জোড়ার দূরত্ব লাগে বলে মেমরি O(n²): কয়েক হাজার সারিতে চলে, লাখে নয়।

৬ · DBSCAN — ঘনত্ব দিয়ে দল, আর noise-কে noise বলা

DBSCAN (Martin Ester, Hans-Peter Kriegel, Jörg Sander আর Xiaowei Xu, KDD ১৯৯৬) দলকে দেখে “ঘন এলাকা” হিসেবে, যাকে ফাঁকা এলাকা আলাদা করে। দুটি parameter:

• eps (ε): প্রতিবেশীর ব্যাসার্ধ।
• min_samples: একটি বিন্দুর ε-এর ভেতরে (নিজেসহ) অন্তত এতগুলো বিন্দু থাকলে সে core point.

পরস্পরের ε-এর মধ্যে থাকা core বিন্দুগুলো জুড়ে একটি দল হয়, অনেকটা BFS দিয়ে graph-এর connected component খোঁজার মতো। কোনো core-এর ε-এর ভেতরে থাকা কিন্তু নিজে core নয় এমন বিন্দু হলো border point, সেই দলে যোগ হয়। আর যে বিন্দু কোনো core-এর কাছেই নেই, সে noise: scikit-learn তাকে label −১ দেয়। k লাগে না, দলের আকৃতি যেকোনো, আর ছিটকে পড়া বিন্দু জোর করে দলে ঢোকে না। দুর্বলতা: দলগুলোর ঘনত্ব আলাদা হলে একটি ε সবার জন্য খাটে না, আর ε ডেটার scale-এর ওপর নির্ভর করে।

নিচের ব্লক দুটি বাঁকা চাঁদের ৩০০টি বিন্দু (M) বানায়, যার আসল দল আমরা জানি; তারপর সঙ্গে ৮টি এলোমেলো ছিটকে পড়া বিন্দু যোগ করে আরেকটি ডেটা (M+8). k-means (আগের ব্লকের নিজের লেখা ফাংশন), AgglomerativeClustering দুই linkage-এ, আর DBSCAN চালিয়ে ARI মাপে (শুধু ৩০০টি চাঁদের বিন্দুতে)। শেষে scipy দিয়ে একই single linkage-এর dendrogram আঁকে।

import matplotlib.pyplot as plt
from sklearn.datasets import make_moons
from sklearn.cluster import AgglomerativeClustering, DBSCAN
from scipy.cluster.hierarchy import linkage, dendrogram, fcluster

M, moon = make_moons(n_samples=300, noise=0.06, random_state=42)    # দুটি বাঁকা চাঁদ
A = np.vstack([M, np.random.default_rng(1).uniform([-1.5, -1], [2.5, 1.5], (8, 2))])  # +৮টি ছিটকে পড়া
def show(name, lab):                               # ARI শুধু ৩০০টি চাঁদের বিন্দুতে
    print(f"{name:17} ARI {adjusted_rand_score(moon, lab[:300]):.3f} | আকার {np.unique(lab, return_counts=True)[1]}")
km = min((kmeans(M, 2, np.random.default_rng(s)) for s in range(10)), key=lambda r: r[2])[0]
show("k-means, M", km)
single = AgglomerativeClustering(n_clusters=2, linkage="single").fit_predict(M)
show("single, M", single)
show("single, M+8", AgglomerativeClustering(n_clusters=2, linkage="single").fit_predict(A))
show("ward, M+8", AgglomerativeClustering(n_clusters=2, linkage="ward").fit_predict(A))
db = DBSCAN(eps=0.2, min_samples=5).fit_predict(A)  # -1 মানে noise
show("DBSCAN, M+8", db)
print("DBSCAN-এর noise:", int((db == -1).sum()), "| তার মধ্যে ছিটকে পড়া:", int((db[300:] == -1).sum()))
L = linkage(M, method="single")                    # scipy: একই single linkage, merge-এর তালিকা
print("scipy বনাম scikit-learn (ARI):", adjusted_rand_score(single, fcluster(L, t=2, criterion="maxclust")))
print("শেষ তিন merge-এর উচ্চতা:", L[-3:, 2].round(3))

fig, axs = plt.subplots(2, 2, figsize=(8, 7))
for ax, P, lab, title in [(axs[0, 0], M, km, "k-means (k=2)"), (axs[0, 1], M, single, "single linkage"),
                          (axs[1, 0], A, db, "DBSCAN on moons + 8 outliers")]:
    ax.scatter(P[:, 0], P[:, 1], c=np.where(lab < 0, 0.5, lab), cmap="coolwarm", s=10)
    ax.set_title(title)
axs[1, 0].scatter(*A[db == -1].T, marker="x", c="k", s=40, label="noise (-1)")
axs[1, 0].legend(loc="lower left")
dendrogram(L, truncate_mode="lastp", p=12, ax=axs[1, 1], color_threshold=0.3, no_labels=True)
axs[1, 1].axhline(0.25, ls="--", c="gray")
axs[1, 1].set_title("single-linkage dendrogram (last 12)")
axs[1, 1].set_ylabel("merge height")
fig.subplots_adjust(hspace=0.35, wspace=0.25)
plt.show()

k-means, M        ARI 0.247 | আকার [151 149]
single, M         ARI 1.000 | আকার [150 150]
single, M+8       ARI 0.000 | আকার [307   1]
ward, M+8         ARI 0.425 | আকার [207 101]
DBSCAN, M+8       ARI 1.000 | আকার [  5 152 151]
DBSCAN-এর noise: 5 | তার মধ্যে ছিটকে পড়া: 5
scipy বনাম scikit-learn (ARI): 1.0
শেষ তিন merge-এর উচ্চতা: [0.124 0.133 0.355]

৭ · PCA-র ধারণা — variance মানে তথ্য, অক্ষ ঘোরানো

এবার দ্বিতীয় কাজ: কলাম কমানো। ধরুন বাজারঘরের প্রতিটি গ্রাহকের ৪০টি কলাম। ছবি আঁকা যায় দুই মাত্রায়, আর অনেক কলামে k-NN ধীর ও দুর্বল হয় (পর্ব ৫৭)। কিন্তু কলামগুলো প্রায়ই একে অন্যের সঙ্গে জড়ানো: যিনি বেশি অর্ডার দেন, তিনি সাধারণত বেশি টাকাও খরচ করেন। তাহলে দুটি কলামের বদলে “কেনাকাটার সক্রিয়তা” নামের একটি কলামেই প্রায় সব কথা বলা যায়।

PCA-র মূল ধারণা: যে দিকে ডেটা বেশি ছড়ানো (variance বেশি), সেদিকে বিন্দুগুলোর পার্থক্য বেশি, তাই সেদিক বেশি তথ্য বহন করে। যে দিকে প্রায় সবাই একই মানে, সেদিক বাদ দিলে কম হারায়। প্রথম principal component (PC1) হলো সেই একক-দৈর্ঘ্যের দিক, যার ওপর বিন্দুগুলোর ছায়ার variance সবচেয়ে বেশি; PC2 হলো PC1-এর সঙ্গে orthogonal দিকগুলোর মধ্যে সবচেয়ে বেশি variance-এর দিক; এভাবে চলে। এই দিকগুলো হলো ডেটার covariance matrix-এর eigenvector, আর প্রতিটির variance তার eigenvalue (পর্ব ১৭-এর “আরও গভীরে”)।

পাঁচটি বিন্দু দিয়ে হাতে, এমন সংখ্যা বেছে যাতে সব হিসাব পরিষ্কার হয়:

৮ · SVD দিয়ে নিজে PCA — explained variance, scaling আর reconstruction

বাস্তবে covariance matrix বানিয়ে eigenvector খোঁজার বদলে সরাসরি centred ডেটার SVD নেওয়া হয় (পর্ব ১৭): X_c = U S V^T. এখানে X_c-এর shape (n, d); V^T-এর প্রতিটি সারি একটি PC-র দিক (d টি ঘর); S-এর কর্ণে singular value s₁ ≥ s₂ ≥ … ; আর প্রতিটি PC-র variance λ_j = s_j² / (n − ১)। কারণ X_c^TX_c / (n − ১) হলো covariance matrix, যা SVD বসালে V (S²/(n − ১)) V^T. SVD পছন্দের কারণ X^TX বানালে নির্ভুলতা কমে (condition number বর্গ হয়, পর্ব ১৭)।

Explained variance ratio = λ_j / ∑λ, অর্থাৎ মোট variance-এর কত ভাগ j-তম PC বহন করে। q টি PC রেখে ফেরত আনলে (X̂ = Z_qV_q^T + μ) আপেক্ষিক reconstruction error = ১ − প্রথম q টি ratio-র যোগফল। নিচের ব্লক অংশ ৩-এর একই RFM ডেটায় (scaled Z আর কাঁচা X) এসব যাচাই করে।

from sklearn.decomposition import PCA
np.set_printoptions(precision=3, suppress=True)

def pca_svd(D):
    Xc = D - D.mean(0)                             # ১. কেন্দ্রে আনা (centering)
    U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
    var = S**2 / (len(D) - 1)                      # প্রতিটি PC-র variance
    return Vt, var / var.sum(), Xc @ Vt.T          # দিক, explained ratio, নতুন স্থানাঙ্ক

Vt, ratio, T = pca_svd(Z)                          # scaled RFM (Z)
sk = PCA().fit(Z)
print("explained variance ratio:", ratio, "| scikit-learn:", sk.explained_variance_ratio_)
print("দিক মেলে (চিহ্ন বাদে)?", np.allclose(np.abs(Vt), np.abs(sk.components_)))
print("PC1-এর ওজন (R, F, M):", Vt[0])
C = np.cov(T.T)                                    # নতুন কলামগুলোর covariance matrix
print("PC-র variance:", C.diagonal(), "| বাকি ঘরের সর্বোচ্চ |মান|:", f"{np.abs(C - np.diag(C.diagonal())).max():.1e}")

Vr, ratio_raw, _ = pca_svd(X)                      # কাঁচা টাকা-দিন-অর্ডার
print("scaling ছাড়া ratio:", [f"{r:.5f}" for r in ratio_raw], "| PC1-এর ওজন:", [f"{v:.4f}" for v in Vr[0]])
print("কাঁচা কলামের σ (R, F, M):", X.std(0).round(1))

for q in (1, 2, 3):                                # q টি PC রেখে ফেরত আনা
    Zhat = T[:, :q] @ Vt[:q] + Z.mean(0)
    err = ((Z - Zhat) ** 2).sum() / ((Z - Z.mean(0)) ** 2).sum()
    print(f"q={q}: reconstruction error {err:.3f} | 1 − cumulative ratio {1 - ratio[:q].sum():.3f}")

explained variance ratio: [0.587 0.283 0.13 ] | scikit-learn: [0.587 0.283 0.13 ]
দিক মেলে (চিহ্ন বাদে)? True
PC1-এর ওজন (R, F, M): [ 0.665 -0.575 -0.476]
PC-র variance: [1.768 0.851 0.391] | বাকি ঘরের সর্বোচ্চ |মান|: 4.5e-16
scaling ছাড়া ratio: ['0.99989', '0.00011', '0.00000'] | PC1-এর ওজন: ['0.0047', '-0.0001', '-1.0000']
কাঁচা কলামের σ (R, F, M): [  93.3    5.5 8045.3]
q=1: reconstruction error 0.413 | 1 − cumulative ratio 0.413
q=2: reconstruction error 0.130 | 1 − cumulative ratio 0.130
q=3: reconstruction error 0.000 | 1 − cumulative ratio 0.000

৯ · load_digits-এ PCA — ছবি, scree plot আর model-এর আগে PCA

scikit-learn-এর সঙ্গে আসা load_digits হাতে লেখা ০ থেকে ৯ অঙ্কের ১৭৯৭টি ছোট ছবি, প্রতিটি ৮ × ৮ পিক্সেল, প্রতিটি পিক্সেলের মান ০ থেকে ১৬; তাই প্রতিটি সারি ৬৪ কলামের (scikit-learn-এর documentation অনুযায়ী এটি UCI Machine Learning Repository-র optical handwritten digits ডেটার একটি অংশ)। প্রশ্ন দুটি: দুই মাত্রায় কেমন দেখায়, আর কম PC-তে classifier-এর কাজ চলে কি না।

Scree plot হলো PC নম্বরের বিপরীতে explained variance (বা তার cumulative যোগফল)-এর রেখা; নামটি পাহাড়ের গোড়ায় জমা পাথরের ঢাল (scree) থেকে, Raymond Cattell-এর ১৯৬৬ সালের কাজে। প্রচলিত নিয়ম: cumulative ratio ৯০% বা ৯৫% ছুঁলে থামা; PCA(n_components=0.95) নিজেই তা করে।

নিচের ব্লকে প্রথমে ডেটা train আর test-এ ভাগ (পর্ব ৫৬), আর PCA শুধু train-এ fit, কারণ test-এর তথ্য দিয়ে দিক বাছলে তা leakage. তারপর q = ২ থেকে ৬৪ PC রেখে make_pipeline(PCA, model) দিয়ে logistic regression (পর্ব ৫৯) আর 5-NN (পর্ব ৫৭)-এর held-out accuracy মাপা; pipeline নিশ্চিত করে PCA প্রতিবার শুধু train-এ fit হয় (বিস্তারিত পর্ব ৬৪, scikit-learn pipeline). 5-NN-এ algorithm="kd_tree" দেওয়া হয়েছে, কারণ ব্রাউজারের Python-এ অন্য algorithm একটি অপ্রাসঙ্গিক warning ছাপে; ফল একই।

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import KNeighborsClassifier
from sklearn.pipeline import make_pipeline

X, y = load_digits(return_X_y=True)                # 1797টি 8×8 হাতের লেখা অঙ্ক, 64 কলাম
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25, stratify=y, random_state=42)
pca = PCA().fit(X_tr)                              # শুধু train-এ fit
cum = np.cumsum(pca.explained_variance_ratio_)
print("shape:", X.shape, "| PC1+PC2:", cum[1].round(3), "| 90% variance-এ লাগে", int(np.argmax(cum >= 0.9)) + 1, "টি PC")

for q in (2, 5, 10, 20, 30, 64):                   # পিক্সেলের একক একই (0–16), তাই scaling নেই
    lr = make_pipeline(PCA(q, random_state=0), LogisticRegression(max_iter=5000)).fit(X_tr, y_tr)
    kn = make_pipeline(PCA(q, random_state=0), KNeighborsClassifier(5, algorithm="kd_tree")).fit(X_tr, y_tr)
    print(f"q={q:2d} | logistic {lr.score(X_te, y_te):.3f} | 5-NN {kn.score(X_te, y_te):.3f}")

P = pca.transform(X_tr)[:, :2]                     # ছবির জন্য train-এর 2-D projection
fig, (a1, a2) = plt.subplots(1, 2, figsize=(9, 4))
sc = a1.scatter(P[:, 0], P[:, 1], c=y_tr, cmap="tab10", s=6)
a1.set_title("load_digits on PC1 and PC2")
a1.set_xlabel("PC1"); a1.set_ylabel("PC2")
fig.colorbar(sc, ax=a1, ticks=range(10), label="digit")
a2.plot(range(1, 65), cum, marker=".")
a2.axhline(0.9, ls="--", c="gray")
a2.set_title("scree: cumulative explained variance")
a2.set_xlabel("number of components"); a2.set_ylabel("cumulative ratio")
fig.subplots_adjust(wspace=0.35, bottom=0.15)
plt.show()

shape: (1797, 64) | PC1+PC2: 0.284 | 90% variance-এ লাগে 21 টি PC
q= 2 | logistic 0.591 | 5-NN 0.620
q= 5 | logistic 0.838 | 5-NN 0.929
q=10 | logistic 0.916 | 5-NN 0.967
q=20 | logistic 0.931 | 5-NN 0.980
q=30 | logistic 0.958 | 5-NN 0.980
q=64 | logistic 0.962 | 5-NN 0.984

এখানে scaling করা হয়নি, কারণ সব পিক্সেলের একক একই (০–১৬ কালির ঘনত্ব)। Scaling করলে কী হয়, তা ডেস্কটপে মাপা হয়েছে:

~/bazarghar-ai$ python digits_scaled.py
q= 5 | scaling ছাড়া 0.838 | StandardScaler-সহ 0.784
q=10 | scaling ছাড়া 0.916 | StandardScaler-সহ 0.858
q=20 | scaling ছাড়া 0.929 | StandardScaler-সহ 0.942
সব সময় 0 এমন পিক্সেল-কলাম: 3

সত্যিকারের run (ডেস্কটপ, Python 3.14.4, scikit-learn 1.9.1; একই ভাগ, PCA → logistic regression). কম PC-তে StandardScaler-সহ accuracy কম (৫টিতে ০.৭৮৪ বনাম ০.৮৩৮), কারণ standardize করলে কিনারার প্রায়-ফাঁকা পিক্সেলও মাঝের পিক্সেলের সমান variance পায়, আর প্রথম PC-গুলো সেই কোলাহলেও টানে। ২০টিতে উল্টো (০.৯৪২ বনাম ০.৯২৯)। নিয়মটি তাই “একক আলাদা হলে scale করুন”; একই এককে দুটোই validation-এ চেষ্টা করুন। (২০ PC-তে ডেস্কটপে ০.৯২৯, ব্রাউজারে ০.৯৩১: version-এর সামান্য সংখ্যাগত পার্থক্য, একটি ছবির।)

t-SNE আর UMAP — শুধু পরিচয়, সতর্কতাসহ

PCA একটি linear পদ্ধতি: নতুন অক্ষ পুরোনো কলামের ওজনযুক্ত যোগ। বাঁকা গঠন (যেমন অংশ ৬-এর চাঁদ) এটি খুলে দেখাতে পারে না। ছবি আঁকার জন্য দুটি nonlinear পদ্ধতি জনপ্রিয়: t-SNE (Laurens van der Maaten আর Geoffrey Hinton, ২০০৮) আর UMAP (Leland McInnes ও সহকর্মীরা, ২০১৮)। দুটিই কাছের বিন্দুকে কাছে রাখতে চায়, আর প্রায়ই PCA-র চেয়ে পরিষ্কার দল দেখায়। কিন্তু:

• দলের মাঝের দূরত্ব অর্থহীন হতে পারে: t-SNE ছবিতে দুটি দল দূরে আছে মানে আসল ডেটায় তারা বেশি আলাদা, এমন বলা যায় না। দলের ছড়ানোও ঘনত্বের কথা বলে না।
• Parameter আর seed-এ ছবি বদলায়: perplexity, n_neighbors বা seed বদলালে ছবি অনেক বদলায়; এলোমেলো ডেটাতেও দল দেখা যেতে পারে।
• Pipeline-এ নয়: scikit-learn-এর TSNE-এ transform নেই, তাই নতুন বিন্দু বসানো যায় না; UMAP আলাদা package (umap-learn).

১০ · AI-এর যুগে — AI assistant-এর লেখা clustering আর PCA কোড যাচাই

“গ্রাহকদের segment করে দিন” বললে AI assistant কোড, ছবি আর দলের সুন্দর নাম দেয়। Unsupervised কাজে সঠিক উত্তর নেই, তাই ভুল কোডও বিশ্বাসযোগ্য দল দেখায়।

AI-এর লেখা clustering বা PCA হাতে পেলে — যাচাইয়ের চেকলিস্ট

• কাঁচা কলামের σ দেখেছি; একক আলাদা হলে scaling হয়েছে, আর তা শুধু train-এ fit (পরে model থাকলে)
• k-means কয়েকটি শুরু থেকে চলেছে (n_init স্পষ্ট); কয়েকটি seed-এ দলের মিল (ARI) মেপেছি
• k বাছাইয়ে elbow আর silhouette দুটোই দেখেছি, আর অন্তত একটি অন্য পদ্ধতির (agglomerative বা DBSCAN) সঙ্গে তুলনা করেছি
• প্রতিটি দলের কয়েকটি সারি নিজে পড়েছি; নাম বর্ণনামূলক আর মানুষকে ছোট করে না
• PCA-র explained ratio আর PC1-এর ওজন দেখেছি; একটি কলাম একাই প্রায় সব নিলে কারণ খুঁজেছি
• কোড নিজে চালিয়েছি, সব warning পড়েছি; কোনো আসল গ্রাহকের ব্যক্তিগত তথ্য AI-তে পেস্ট করিনি