قراءة معمقة في SVD

قمة الجبل

قد يبدو الحديث عن SVD في عالم الجبر الخطي كالقفز مباشرةً إلى قمة جبل. لماذا نبدأ من أعلى النظرية قبل أن نعرف الطريق إليها؟

لا تقلق. لن نحفظ القمة ولن نضيع وسط متاهات البراهين. سنفعل شيئًا أبسط: سنقف فوقها قليلًا، ثم نبدأ في تفكيكها حتى نرى ممَّا تكوَّن الجبل.

من التجميع إلى التفكيك

وبما أننا نتحدث عن التفكيك — Decomposition، فلنبدأ من الاتجاه المعاكس: التجميع — Composition.

أمامنا ثلاث مصفوفات، $A$ و$B$ و$C$. أبعادها متوافقة، ولذلك نستطيع أن نصلها ببعضها:

هذه هي الفكرة في أبسط صورة: نعرف القطع، نجمعها، فيظهر شيء جديد. وإذا دخل متجه $x$ إلى هذه السلسلة فإن العمل يبدأ من اليمين: $Wx=ABCx=A(B(Cx))$، أي

لنأخذ مثالًا صغيرًا:

وعند جمع هذه التحويلات نحصل على

التجميع في Python

1import numpy as np2 3A = np.array([[1., 2.],4              [0., 1.]])5 6B = np.array([[2., 0.],7              [0., 0.5]])8 9C = np.array([[0., -1.],10              [1.,  0.]])11 12W = A @ B @ C13print(W)

هذا هو الاتجاه السهل في القصة: القطع معروفة، ولذلك نبني منها المصفوفة النهائية.

هنا تلمس التجميع بدل أن تقرأ تعريفه فقط. غيّر $x$، واضغط على المصفوفات، وراقب كيف تمر المتجهة من $C$ إلى $B$ ثم إلى $A$.

والآن اقلب السؤال: بدل أن أعطيك القطع، أضع أمامك مصفوفة واحدة اسمها $D$. نرى صفوفها وأعمدتها، ونستطيع أن نحسب رتبتها وكل ما يظهر لنا منها. ثم أقول لك: افتحها. أرني البنية التي تختبئ في داخلها.

لو لم تكن لدينا أدوات، لاستطعنا التخمين حتى تحترق النجوم. لكن هذا بالضبط هو السبب الذي جعل علماء الجبر الخطي يبنون نظريات الـdecomposition: لا نريد أي عوامل عشوائية؛ نريد قطعًا لها معنى وبنية نستطيع الاستفادة منها.

قد تكون القطع مثلثية كما في LU، أو اتجاهات orthonormal ومعاملات مثلثية كما في QR، أو عاملًا مثلثيًا وانعكاسه كما في Cholesky. وقد نغيّر زاوية النظر نفسها كما في Spectral decomposition وSchur، أو نفصل الدوران عن التمدد كما في Polar decomposition. ثم نصل إلى SVD، الذي سيأخذ المسرح من هنا.

أطلس التفكيك

الأسماء كثيرة، لكن الفكرة واحدة: لا تقاتل المصفوفة بالشكل الذي وصلت به إليك. ابحث عن تمثيل يجعل عمل كل قطعة واضحًا.

في العنصر التالي اختر LU أو QR أو Cholesky أو Spectral أو SVD. افتح العوامل، ثم أعد بناء المصفوفة الأصلية. لا تحفظ الصيغ؛ راقب النمط وهو يتكرر.

أطلس تفكيك صغير في Python

1import numpy as np2 3A = np.array([[4., 2.],4              [2., 3.]])5 6# QR7Q, R = np.linalg.qr(A)8print(np.allclose(A, Q @ R))9 10# SVD11U, s, Vt = np.linalg.svd(A)12Sigma = np.diag(s)13print(np.allclose(A, U @ Sigma @ Vt))14 15# Spectral decomposition for symmetric A16lam, E = np.linalg.eigh(A)17print(np.allclose(A, E @ np.diag(lam) @ E.T))18 19# Cholesky for positive-definite A20L = np.linalg.cholesky(A)21print(np.allclose(A, L @ L.T))

وباستخدام SciPy يمكن حساب بعض التفكيكات مباشرة:

1import numpy as np2from scipy.linalg import lu, schur, polar3 4A = np.array([[4., 2.],5              [2., 3.]])6 7P, L, U_lu = lu(A)8T, Z = schur(A)9Qp, H = polar(A)10 11print(np.allclose(A, P @ L @ U_lu))12print(np.allclose(A, Z @ T @ Z.T))13print(np.allclose(A, Qp @ H))

ومن هنا لا يعود SVD قفزة غريبة. أصبح مجرد التفكيك الذي سنقترب منه أكثر من أي تفكيك آخر.

الجوهر فقط

لأي مصفوفة حقيقية $A\in\mathbb R^{m\times n}$، يمكننا كتابة

في الصورة الكاملة يكون $U$ مصفوفة متعامدة من الحجم $m\times m$، و$V$ مصفوفة متعامدة من الحجم $n\times n$، بينما $\Sigma$ مصفوفة قطرية مستطيلة من الحجم $m\times n$.

على قطر $\Sigma$ توجد القيم المفردة:

أعمدة $V$ هي right singular vectors، وأعمدة $U$ هي left singular vectors. ولو كانت العناصر مركبة نستبدل $T$ بالـconjugate transpose $^*$.

هذه هي الصيغة. أمّا المعنى فيظهر عندما نرى وظيفة كل قطعة.

احسب SVD وتحقق منه

1import numpy as np2 3A = np.array([[3., 1.],4              [1., 2.]])5 6U, s, Vt = np.linalg.svd(A)7Sigma = np.diag(s)8 9print("U =\n", U)10print("singular values =", s)11print("V^T =\n", Vt)12 13print("U^T U =\n", U.T @ U)14print("V^T V =\n", Vt @ Vt.T)15 16A_reconstructed = U @ Sigma @ Vt17print(A_reconstructed)18print("error =", np.linalg.norm(A - A_reconstructed))

ولمصفوفة مستطيلة استخدم

1import numpy as np2 3A = np.array([[3., 1.],4              [1., 2.]])5U, s, Vt = np.linalg.svd(A, full_matrices=False)

إذا كانت نسخة compact/economy هي كل ما تحتاجه.

من اليمين تبدأ الحكاية

إذا أدخلنا متجهة $x$، فلدينا $Ax=U\Sigma V^Tx$.

ابدأ من $V^T$. هذه القطعة لا “تشوّه” المتجهة كيفما اتفق؛ هي تسأل: كم تحتوي $x$ من كل اتجاه خاص $v_i$؟ فعلًا، إذا كتبنا $V=[v_1\ v_2\ \cdots\ v_n]$، فإن

إذن $V^T$ يحلل المدخل في إحداثيات صُممت خصيصًا لهذه المصفوفة.

ثم تصل الإحداثيات إلى $\Sigma$. لا توجد هنا مزاوجة معقدة بين المحاور؛ كل اتجاه يُضرب في رقم واحد:

وهنا تكمن بساطة SVD: التحويل المعقد صار في هذه الإحداثيات مجرد stretch/shrink مستقل لكل محور.

وأخيرًا تأتي $U$ لتضع هذه المركبات بعد تمددها في اتجاهاتها المناسبة داخل فضاء الخرج.

لذلك يمكن أن نحفظ القصة لا الحروف:

يحلّل (V^T) المدخل، ثم تمدّد (Sigma) مركباته أو تقلّصها، وأخيرًا يوجّه (U) الخرج.

والجملة التي تختصر كل ذلك هي

أدخل الاتجاه الخاص $v_i$، فتخرج في اتجاه $u_i$ بعد أن يتغير طولها بمقدار $\sigma_i$.

سر الجمال

هنا يظهر أول سبب يجعل SVD جميلًا: الـeigendecomposition التقليدي يريد مصفوفة مربعة، لأن معادلة

تحتاج أن يعيش $v$ و$Av$ في الفضاء نفسه.

لكن إذا كانت $A:\mathbb R^n\to\mathbb R^m$ و$m\neq n$، فليس هناك سبب لأن يكون فضاء المدخل هو فضاء الخرج.

SVD لا يحاول إجبار الفضاءين على أن يكونا واحدًا. يعطي المدخل أساسه الخاص $V$، ويعطي الخرج أساسه الخاص $U$. هذا ليس حلًا التفافيًا للمشكلة؛ هذه هي الفكرة التي تجعل SVD طبيعيًا للمصفوفات المستطيلة.

حين تتحول الدائرة

تخيل دائرة الوحدة في بعدين. كل نقطة عليها تمثل متجهة طولها واحد. عندما نطبّق $V^T$، تتغير الإحداثيات دون أن تتغير الأطوال لأن $V$ orthogonal. ثم تأتي $\Sigma$ فتمد محورًا أكثر من الآخر، فتحول الدائرة إلى ellipse. بعدها تدور $U$ النتيجة أو تعكسها لتضعها في الاتجاه النهائي.

دع الهندسة تتحرك

1import numpy as np2import matplotlib.pyplot as plt3 4A = np.array([[3., 1.],5              [1., 2.]])6 7U, s, Vt = np.linalg.svd(A)8Sigma = np.diag(s)9 10theta = np.linspace(0, 2*np.pi, 500)11circle = np.vstack([np.cos(theta), np.sin(theta)])12 13after_vt = Vt @ circle14after_sigma = Sigma @ after_vt15after_u = U @ after_sigma16 17plt.figure(figsize=(6, 6))18plt.plot(after_u[0], after_u[1])19plt.axhline(0)20plt.axvline(0)21plt.gca().set_aspect("equal", adjustable="box")22plt.title("Final image of the unit circle under A")23plt.show()

للتعليم من الأفضل متابعة circle ثم after_vt ثم after_sigma ثم after_u مرحلةً بعد مرحلة. وفي النسخة التفاعلية نتابع الجسم نفسه وهو يتغير بدل الاعتماد على صور منفصلة.

لهذا يقال كثيرًا إن SVD هو “دوران/انعكاس، ثم تمدد، ثم دوران/انعكاس”. العبارة مفيدة، لكن النسخة الأدق هي أن $V^T$ يختار نظام إحداثيات المدخل، و$\Sigma$ يقوم بالفعل المستقل على كل محور، و$U$ يختار نظام إحداثيات الخرج.

من أين تأتي هذه القطع؟

الآن وصلنا إلى سؤال مهم: هل اخترعنا هذه الاتجاهات لأنها تبدو جميلة؟ لا. ابدأ من المعادلة $Av_i=\sigma_i u_i$.

اضرب من اليسار بـ$A^T$:

وفي SVD نحصل أيضًا على $A^Tu_i=\sigma_i v_i$، إذن $A^TAv_i=\sigma_i^2v_i$.

ها هو السر: الـright singular vectors هي eigenvectors للمصفوفة $A^TA$، ومربعات singular values هي eigenvalues لها.

وبالمثل $AA^Tu_i=\sigma_i^2u_i$، فتكون left singular vectors eigenvectors لـ$AA^T$.

وهذا يفسر لماذا القيم المفردة غير سالبة: $A^TA$ موجبة شبه محددة، وبالتالي eigenvalues الخاصة بها غير سالبة، ثم نأخذ الجذر $\sigma_i=\sqrt{\lambda_i}$.

أول خيط: $A^TA$

إذا أردت أن تبني SVD للتعلم لا للحساب العددي الاحترافي، ابدأ بـ$A^TA$. استخرج eigenvectors متعامدة معيارية $v_i$ ورتب eigenvalues تنازليًا. من كل eigenvalue $\lambda_i$ اصنع

وحين تكون $\sigma_i\neq0$، مرر $v_i$ داخل $A$، ثم اقسم على مقدار التمدد:

بهذه الطريقة لا نسحب $U$ و$\Sigma$ و$V$ من قبعة ساحر. $A^TA$ تخبرنا بالاتجاهات المهمة في فضاء المدخل، eigenvalues تخبرنا بقوة التمدد، و$A$ نفسها تخبرنا أين تصل هذه الاتجاهات في فضاء الخرج.

هذا الطريق ممتاز للفهم. أما في العمل العددي الحقيقي فسنعتمد على خوارزميات SVD مستقرة مباشرة بدل تكوين $A^TA$ دائمًا، لأن تربيع condition number قد يزيد المشكلات العددية.

ابنِ SVD من $A^TA$

هذا مفيد لفهم النظرية. أما في العمل العددي الحقيقي فالأفضل استخدام np.linalg.svd مباشرة.

1import numpy as np2 3def svd_from_ata(A, tol=1e-12):4    A = np.asarray(A, dtype=float)5 6    # 1. Symmetric positive-semidefinite Gram matrix7    G = A.T @ A8 9    # 2. Its orthonormal eigenvectors become right singular vectors10    eigenvalues, V = np.linalg.eigh(G)11 12    # 3. Sort from largest to smallest13    order = np.argsort(eigenvalues)[::-1]14    eigenvalues = np.maximum(eigenvalues[order], 0.0)15    V = V[:, order]16 17    # 4. Singular values are square roots18    s = np.sqrt(eigenvalues)19 20    # 5. Keep the nonzero part for the condensed SVD21    keep = s > tol22    s = s[keep]23    V = V[:, keep]24 25    # 6. Av_i = sigma_i u_i26    U = A @ V / s27 28    return U, s, V.T

برهان الوجود

ما الذي نريد إثباته فعلًا؟

حتى الآن عرفنا كيف نفهم SVD وكيف نبنيه بطريقة تعليمية. لكن تبقى قفزة مهمة: لماذا نضمن أن كل مصفوفة تملك SVD أصلًا؟

لنفترض $A\in\mathbb C^{m\times n}$، و$r=\operatorname{rank}(A)$.

نريد أن نثبت وجود صيغة condensed SVD:

بحيث أعمدة $X$ وأعمدة $Y$ orthonormal، و

إذا كانت الرموز المركبة مزعجة، تذكر فقط أن $A^*=A^T$ في الحالة الحقيقية.

الأداة التي سنستعملها هي spectral theorem: كل مصفوفة Hermitian تملك أساسًا orthonormal من eigenvectors، ويمكن كتابتها في صورة

المشكلة الوحيدة أن $A$ قد لا تكون مربعة ولا Hermitian. إذن بدل أن نحاول إجبار $A$ على شيء ليست عليه، سنبني مصفوفة أكبر تحمل $A$ داخلها.

الحيلة التي تفتح الباب

عرّف

حجمها $(m+n)\times(m+n)$، فهي مربعة. والأهم أن $W^*=W$، إذن $W$ Hermitian، وبالتالي spectral theorem مفتوح لنا.

خذ eigenvector لـ$W$ واكتبه على قطعتين:

بتوسيع الضرب الكتلي نحصل على

أي $Ay=\sigma x$ و$A^*x=\sigma y$.

انظر إلى أول معادلة. إنها تقريبًا التعريف الذي كنا نبحث عنه: $y$ يتصرف كـright singular vector، و$x$ كـleft singular vector، و$\sigma$ كـsingular value.

افحص برهان الـblock matrix عدديًا

1import numpy as np2 3A = np.array([[3., 1.],4              [0., 2.],5              [2., 2.]])6 7U, s, Vt = np.linalg.svd(A, full_matrices=False)8 9m, n = A.shape10W = np.block([11    [np.zeros((m, m)), A],12    [A.T, np.zeros((n, n))]13])14 15# Hermitian/symmetric, so eigvalsh is appropriate16w_eigenvalues = np.linalg.eigvalsh(W)17 18expected = np.sort(19    np.concatenate([20        -s,21        np.zeros(m + n - 2*len(s)),22        s23    ])24)25 26print(w_eigenvalues)27print(expected)28print(np.allclose(w_eigenvalues, expected))

هذه التجربة العددية تعرض ما يتنبأ به البرهان تمامًا:

هذه هي الحيلة المركزية للبرهان كله: نحول مشكلة SVD لمصفوفة ربما تكون مستطيلة إلى eigenvalue problem لمصفوفة Hermitian أكبر.

لماذا تأتي القيم أزواجًا؟

إذا كان

eigenvector بقيمة $+\sigma$، فانظر إلى

باستخدام $Ay=\sigma x$ و$A^*x=\sigma y$ نجد

إذن القيم غير الصفرية تأتي في أزواج $+\sigma_i,-\sigma_i$.

وبما أن $W$ Hermitian، فإن eigenvectors المرتبطة بقيم eigenvalues مختلفة يمكن اختيارها متعامدة. هذه المعلومة هي التي ستعطينا لاحقًا orthonormality للاتجاهات المفردة.

نجمع القطع

نختار

بحيث يكون $z_i^*z_i=2$.

قد يبدو الرقم 2 غريبًا، لكنه مقصود: نريد في النهاية أن يكون لكل من $x_i$ و$y_i$ طول يساوي 1.

لدينا $x_i^*x_i+y_i^*y_i=2$.

ومن تعامد eigenvector المرتبط بـ$+\sigma_i$ مع النظير المرتبط بـ$-\sigma_i$ نحصل على

نجمع المعادلتين فنحصل على $x_i^*x_i=1$، ونطرحهما فنحصل على $y_i^*y_i=1$.

الآن نجمع المتجهات في مصفوفتين:

ونضع القيم المفردة في $\Sigma_r=\operatorname{diag}(\sigma_1,\ldots,\sigma_r)$.

وباستخدام spectral decomposition للمصفوفة $W$، ثم ضرب الكتل ومقارنة الجزء العلوي الأيمن، نحصل على ما أردناه:

لماذا تبقى الاتجاهات نظيفة؟

وحدة الطول حصلنا عليها بالفعل. بقي أن نثبت أن الأعمدة المختلفة متعامدة.

لـ$i\neq j$، تعامد eigenvectors لـ$W$ يعطينا

وعندما نقارن eigenvector الموجب للأول مع النظير ذي الإشارة السالبة للثاني نحصل على

اجمع المعادلتين: $2x_i^*x_j=0 \Longrightarrow x_i^*x_j=0$. واطرحهما: $2y_i^*y_j=0 \Longrightarrow y_i^*y_j=0$. وهكذا $X^*X=I_r$ و$Y^*Y=I_r$.

لم نعد نملك حدسًا فقط. لقد أثبتنا وجود الـcondensed SVD.

البرهان في نفس واحد

إذا ضاعت التفاصيل، لا تحفظ الصفحات. احفظ الحركة.

نبدأ بـ$A$، ثم نخفيها داخل

$W$ Hermitian، لذلك تعطينا spectral theorem eigenvectors متعامدة معيارية. نقسم كل eigenvector إلى جزأين $x_i$ و$y_i$، فتخرج لنا المعادلتان

نجمع $x_i$ في $X$، و$y_i$ في $Y$، والقيم $\sigma_i$ في $\Sigma_r$، فنستعيد

الـsingular vectors لم تظهر من العدم. كانت مختبئة داخل eigenvectors لمصفوفة Hermitian صممناها بعناية.

من حاصل ضرب إلى طبقات

حتى الآن رأينا SVD كآلة من ثلاث مراحل. الآن سنديرها قليلًا لنرى صورة ثانية لا تقل أهمية. اكتب $U=[u_1\ u_2\ \cdots]$ و$V=[v_1\ v_2\ \cdots]$.

لأن $\Sigma$ قطرية، يمكن فتح حاصل الضرب إلى مجموع:

كل $u_iv_i^T$ مصفوفة rank-1، لأن كل أعمدتها مضاعفات عددية من $u_i$.

إذن SVD يقول شيئًا مدهشًا وبسيطًا في الوقت نفسه:

المصفوفة طبقات rank-1 مرتبة من الأقوى إلى الأضعف.

انظر إلى المصفوفة

حتى heatmap بسيطة للمصفوفة قد تكفي لكي تصبح البنية مرئية:

1import numpy as np2import matplotlib.pyplot as plt3 4A = np.array([5    [5., 4., 3., 2., 1.],6    [4., 3.2, 2.4, 1.6, 0.8],7    [3., 2.4, 1.8, 1.2, 0.6],8    [2., 1.6, 1.2, 0.8, 0.4],9    [1., 0.8, 0.6, 0.4, 0.2],10])11 12plt.figure(figsize=(7, 5))13im = plt.imshow(A, aspect="auto")14plt.colorbar(im)15plt.title("Matrix values")16plt.xlabel("column")17plt.ylabel("row")18plt.show()

استعمل العرض نفسه كسرد متصل: ابدأ بالمصفوفة الأصلية $A$، ثم اعزل طبقة rank-1، ثم اعرض truncated reconstruction $A_k$، وفي النهاية انظر إلى residual $A-A_k$. هكذا يصبح matrix decomposition شيئًا نراه يتحرك بدل أن يكون مجرد رموز.

اسحب طبقات rank-1 واحدةً واحدة

1import numpy as np2 3A = np.array([4    [5., 4., 3., 2., 1.],5    [4., 3.2, 2.4, 1.6, 0.8],6    [3., 2.4, 1.8, 1.2, 0.6],7    [2., 1.6, 1.2, 0.8, 0.4],8    [1., 0.8, 0.6, 0.4, 0.2],9])10U, s, Vt = np.linalg.svd(A, full_matrices=False)11layers = [s[i] * np.outer(U[:, i], Vt[i, :]) for i in range(len(s))]12 13# Exact reconstruction14A_again = sum(layers)15print(np.allclose(A, A_again))

كل layer_i رتبتها لا تتجاوز 1.

ولعرض إحدى الطبقات:

1import matplotlib.pyplot as plt2import numpy as np3 4A = np.array([5    [5., 4., 3., 2., 1.],6    [4., 3.2, 2.4, 1.6, 0.8],7    [3., 2.4, 1.8, 1.2, 0.6],8    [2., 1.6, 1.2, 0.8, 0.4],9    [1., 0.8, 0.6, 0.4, 0.2],10])11U, s, Vt = np.linalg.svd(A, full_matrices=False)12first_layer = s[0] * np.outer(U[:, 0], Vt[0, :])13 14plt.figure()15plt.imshow(first_layer, aspect="auto")16plt.colorbar()17plt.title("First rank-1 SVD layer")18plt.show()

القيمة $\sigma_i$ هي قوة الطبقة. إذا كانت صفرًا، فالطبقة لا تضيف شيئًا. ومن هنا نحصل على حقيقة مهمة:

كم قطعة نحتاج؟

يمكن أيضًا فهم rank من زاوية البناء. إذا كانت

فيمكن بناء $A$ من $k$ قطع rank-1، ولا يمكن فعل ذلك باستخدام $k-1$ قطع rank-1 فقط.

هذه ليست مجرد صورة جميلة؛ إنها تعيد تفسير معنى الرتبة نفسه: كم اتجاهًا مستقلًا نحتاج حتى نبني المصفوفة؟

وهذا يقود طبيعيًا إلى rank factorization. إذا

وكان لـ$Y$ عدد $k$ من الأعمدة، فإن كل أعمدة $A$ تعيش داخل span لتلك الأعمدة، ولذلك

إذا كانت rank الفعلية تساوي $r$، فلن تستطيع ضغط البعد الداخلي تحت $r$ مع الحفاظ على المساواة التامة.

SVD يعطينا rank factorization خاصة جدًا، لكنها لا تكتفي بأي أساس؛ تختار اتجاهات orthonormal وترتبها حسب القوة.

حين نقبل أن نفقد قليلًا

لنفترض أن

إذا كانت أول عدة قيم مفردة كبيرة والبقية صغيرة، قد نقرر الاحتفاظ بأول $k$ طبقات فقط:

نحن هنا لا نزعم أن $A_k=A$. نحن نختار أن نخسر جزءًا من المعلومة مقابل تمثيل أبسط.

truncated SVD

1import numpy as np2 3U, s, Vt = np.linalg.svd(A, full_matrices=False)4 5k = 36A_k = (U[:, :k] * s[:k]) @ Vt[:k, :]7 8error = np.linalg.norm(A - A_k, ord="fro")9tail_error = np.sqrt(np.sum(s[k:]**2))10 11print("Frobenius error:", error)12print("tail singular-value formula:", tail_error)

القيمتان تتفقان حتى حدود الدقة العددية.

ما الذي خسرناه؟ الفرق هو $A-A_k$.

ولقياس حجمه نستعمل Frobenius norm:

وبسبب تعامد طبقات SVD نحصل على

أي أن خطأ الحذف ليس لغزًا؛ إنه مجموع طاقات الطبقات التي قررنا تركها خلفنا.

لماذا القطع الأولى هي الأفضل؟

هنا تأتي واحدة من أجمل نتائج القصة. truncated SVD ليس مجرد طريقة معقولة لاختيار rank-$k$ approximation. إنه الأفضل في Frobenius norm.

إذا كان $B$ أي مصفوفة تحقق $\operatorname{rank}(B)\le k$، فإن

أي أنه إذا أعطيتك ميزانية لا تسمح إلا بـ$k$ اتجاهات مستقلة، فلن تجد rank-$k$ matrix أقرب إلى $A$ من truncated SVD.

الفكرة وراء البرهان هي أن row space لأي $B$ رتبتها لا تتجاوز $k$ لا يمكن أن تلتقط أكثر من $k$ اتجاهات مستقلة. إذا أسقطنا $A$ على أي فضاء بعده $k$، فإن أكبر طاقة يمكن الاحتفاظ بها تأتي من الاتجاهات المرتبطة بأكبر singular values. لذلك اختيار $v_1,\ldots,v_k$ ليس ذوقًا؛ هو الاختيار الذي يحتفظ بأكبر مقدار ممكن من $\sum\sigma_i^2$.

وبالتالي يكون الجزء المفقود أصغر ما يمكن:

هذه هي روح Eckart–Young: إذا كان عليك أن تنسى، فانْسَ أضعف الاتجاهات أولًا.

أين نتوقف؟

النظر إلى singular-value spectrum يشبه النظر إلى طبقات الجبل من جانبه. إذا هبطت القيم بسرعة ثم ظهر ذيل طويل صغير، فهذا يخبرنا أن كثيرًا من البنية يمكن وصفه بعدد قليل من الاتجاهات.

singular values وخطأ التقريب

1import numpy as np2 3A = np.array([[3., 1.],4              [0., 2.],5              [2., 2.]])6U, s, Vt = np.linalg.svd(A, full_matrices=False)7 8k = 39A_k = (U[:, :k] * s[:k]) @ Vt[:k, :]10 11error = np.linalg.norm(A - A_k, ord="fro")12tail_error = np.sqrt(np.sum(s[k:]**2))13 14print("Frobenius error:", error)15print("tail singular-value formula:", tail_error)

خطأ Frobenius الأمثل عند rank-$k$:

1import numpy as np2import matplotlib.pyplot as plt3 4A = np.array([[3., 1.],5              [0., 2.],6              [2., 2.]])7s = np.linalg.svd(A, compute_uv=False)8 9errors = []10for k in range(len(s) + 1):11    errors.append(np.sqrt(np.sum(s[k:]**2)))12 13plt.figure()14plt.plot(range(len(s)+1), errors, marker="o")15plt.xlabel("k")16plt.ylabel("||A - A_k||_F")17plt.title("Best rank-k approximation error")18plt.show()

مقياس شائع هو الطاقة المحتفظ بها:

لا يوجد $k$ سحري يصلح لكل مشكلة. اختيار $k$ هو مقايضة بين البساطة والدقة، وفي التطبيقات الحقيقية قد تدخل أيضًا قيود الذاكرة والزمن وأداء المهمة نفسها.

وهنا تدخل PCA

لنفترض أن $X$ data matrix مركز حول المتوسط، وفيه $n$ observations على الصفوف و$p$ features على الأعمدة.

نفككه:

مصفوفة covariance هي

باستخدام SVD:

إذن

وهنا نرى الجسر بوضوح: أعمدة $V$ هي principal directions، وeigenvalues للـcovariance تساوي

إذن SVD لا يصبح PCA تلقائيًا لمجرد أننا نفكك أي مصفوفة. نحتاج أولًا إلى تفسير data matrix، وعادةً إلى mean centering، ثم نستخدم SVD كآلة حسابية ونظرية لاستخراج اتجاهات أكبر variance.

PCA مباشرة من SVD

1import numpy as np2 3X = np.array([4    [2.5, 2.4], [0.5, 0.7], [2.2, 2.9], [1.9, 2.2], [3.1, 3.0],5    [2.3, 2.7], [2.0, 1.6], [1.0, 1.1], [1.5, 1.6], [1.1, 0.9],6])7 8# X: rows = observations, columns = features9X = np.asarray(X, dtype=float)10 11# 1. Center each feature12mean = X.mean(axis=0, keepdims=True)13Xc = X - mean14 15# 2. SVD of centered data16U, s, Vt = np.linalg.svd(Xc, full_matrices=False)17 18# Principal directions19components = Vt20 21# Principal-component scores22scores = Xc @ components.T23# equivalently: scores == U * s24 25# Covariance eigenvalues / explained variances26n = X.shape[0]27explained_variance = s**2 / (n - 1)28explained_ratio = explained_variance / explained_variance.sum()29 30print("components:\n", components)31print("explained variance:", explained_variance)32print("explained ratio:", explained_ratio)

ولخفض البيانات إلى $k$ components: