ניתוח אשכולות ספקטרלי

ניתוח אשכולות ספקטרלי הוא שיטה לפתרון בעיית ניתוח אשכולות באמצעות שימוש בערכים העצמיים (ספקטרום) של מטריצת הלפלסיאן המתארת את המידע. במקרה ובו הקלט ניתן כגרף, מטריצת הלפלסיאן נחשבת כקלט של האלגוריתם (מאחר שקיימת התאמה אחת ויחידה בין גרף בעל $n$ צמתים לבין מטריצת לפלסיאן מסדר $n\times n$ ), אמנם במקרים רבים הקלט אינו מתקבל כך, (אלא למשל אוסף נקודות במימד $R^{m}$ ) ונדרשת בניית גרף כשלב מקדים לשימוש באלגוריתם. קודקודי הגרף מייצגים את המידע אותו רוצים לנתח, וקשתות הגרף מייצגות את הקשר בין הפרטים בגרף.

הגדרת הבעיה

תיאור הבעיה

ניתוח המידע לאשכולות שקול לחלוקת הגרף לתת קבוצות זרות זו לזו, נרצה למצוא חלוקה שאכן מחלקת את המידע בצורה "טבעית" כלשהי, כלומר אינטואיטיבית קל להבין כי נרצה למצוא חלוקה כזו שבה הקשרים בין חברי קבוצות שונות רופפים עד לא קיימים, ובנוסף הקשרים שבין חברי כל קבוצה זה לזה חזקים, דרך נוספת שעל פיה ניתן לחלק את המידע היא אוסף קבוצות אשר הקשרים ביניהן רופפים, וכן גודלן דומה זה לזה.

יש לשים לב כי בעוד התנאי הראשון הוא בעצם בעיית החתך המינימלי (אנ'), אשר ידועה כבעיה הניתנת לפתרון בזמן ריצה פולינומי, קיומו אינו מספיק לפתרון אופטימלי לבעיית הניתוח מאחר שפעמים רבות תוצאתו היא פיצול קודקוד אחד (או קבוצה קטנה של קודקודים) משאר הגרף, ובהוספת התנאי השני מתקבלת בעיה NP-קשה^[1]. קירוב לפתרון בעיה זו ניתן להיעשות על ידי שימוש במטריצת הלפלסיאן של הגרף.

הגדרה פורמלית

בצורה פורמלית ניתן לתאר את הבעיה באחת מהדרכים הבאות - בהינתן גרף ממושקל ולא מכוון $G=(V,E,w)$ , ומספר $k$ נרצה למצוא חלוקה המביאה למינימום את אחת משתי פונקציות המטרה הבאות:

$RatioCut(A_{1},...,A_{k}):=\sum _{i\,=\,1}^{k}{\frac {cut(A_{i},{\overline {A_{i}}})}{|A_{i}|}}$

$nCut(A_{1},...,A_{k}):=\sum _{i\,=\,1}^{k}{\frac {cut(A_{i},{\overline {A_{i}}})}{vol(A_{i})}}$

כאשר מגדירים:

$A_{i}$ היא תת-קבוצה של צמתים בגרף, מתקיים $A_{i}\cap A_{j}=\emptyset$ , $\cup _{i=1}^{n}A_{i}=V$
$|A_{i}|$ הוא מספר הצמתים ברכיב $A_{i}$ .
$vol(A_{i})$ הוא סכום ערכי הקשתות המחברות את צמתי $A_{i}$ זה לזה.
$cut(A_{i},{\overline {A_{i}}})$ הוא סכום ערכי הקשתות המחברות בין הקבוצה $A_{i}$ לבין שאר הגרף.
$k$ הוא מספר האשכולות, כלומר מספר תת-הקבוצות שאליהן נחלק את הגרף.

ההבדל בין שתי הפונקציות הוא שבעוד $RatioCut$ "מענישה" רכיבי קשירות בעלי מספר צמתים נמוך, ובשימוש בה נצפה לתוצאה של אשכולות בעלי מספר צמתים דומה, $nCut$ "מענישה" רכיבי קשירות בעלי סכום קשתות נמוך, ולכן נצפה לקשרים חזקים בתוך כל אשכול.

בדוגמה שבתמונה ניתן להראות על ידי חישוב כי שתי פונקציות המטרה יביאו לחלוקה שונה, שימוש בפונקציה $RatioCut$ יוביל ל"שבירת" הקשת שערכה 14, לצורך חלוקת הגרף לשתי הקבוצות (השוות בגודלן) $\{1,2,3\},\{4,5,6\}$ . לעומת זאת, השימוש בפונקציה $nCut$ יוביל לחלוקה לקבוצות $\{1,2\},\{3,4,5,6\}$ , מאחר שגודל הקבוצות לא רלוונטי אלא סכום הקשתות, ותחת חישוב זה שבירת הקשת בעלת משקל 13 עדיפה על הקשת הכבדה יותר.

ניתן להניח כי ההבדלים בין הגישות נובעים מסיבות היסטוריות, $RatioCut$ פורסמה בהקשר של חלוקה של רכיבי מעגל חשמלי לתת קבוצות^[2] אשר נוח "לטפל" בכל אחת בפני עצמה, בהקשר זה חשוב לדאוג שהקבוצות תהיינה בעלות גודל דומה וכן שלא יהיו קשרים חזקים בין קבוצה לקבוצה (לדוגמה, כדי שיהיה אפשר למקם כל תת-קבוצה במקום בלתי תלוי ככל האפשר בשאר המעגל החשמלי). לעומתה, $nCut$ פורסמה בהקשר של סגמנטציה^[3], ובהקשר זה חשוב לשים לב לכלול באותו אשכול אך ורק פרטים הקשורים זה לזה.

מטריצת הלפלסיאן

הגדרת מטריצת הלפלסיאן

מטריצת הלפלסיאן מוגדרת בצורה הבאה:

$L:=D-A$

כאשר $A$ היא מטריצה השכנות של הגרף (במקרה של גרף לא ממושקל) או מטריצת הדמיון(אנ') של הגרף (במקרה של גרף ממושקל), וכן $D$ היא מטריצה אלכסונית שבה האיבר $D_{ii}$ הוא הדרגה של האיבר ה $i$ , כלומר מתקיים $\sum _{i\,=\,1}^{n}A_{ij}=D_{jj}$

למטריצת הלפלסיאן קיימות שתי גרסאות מנורמלות:

מטריצת לפלסיאן סימטרית מנורמלת (symmetric normalized Laplacian matrix) מוגדרת כך:

$L^{\text{sym}}:=D^{-1/2}LD^{-1/2}=I-D^{-1/2}AD^{-1/2}$

מטריצת לפלסיאן מנורמלת של מהלך אקראי (random-walk normalized Laplacian matrix) מוגדרת כך: $L^{\text{rw}}:=D^{-1}L=I-D^{-1}A$ תוצאת הנרמול הי מטריצה שבה ערכי האלכסון הם 1, והשימוש במטריצה מנורמלת ייעשה כדי להביא למינימום את פונקציית המטרה $ncut$ , יש לציין כי בגרף רגולרי בדרך כלל לא נדרש נרמול, שכן תוצאת השימוש מטריצה הלא מנורמלת יהיה דומה מאוד למנורמלת.

שימוש במטריצת הלפלסיאן

מתכונות מטריצת הלפלסיאן הלא מנורמלת, ניתן להוכיח^[4] כי בהינתן תת-הגרף $A_{j}$ וכן וקטור אינדיקטור(אנ') מנורמל $h$ של תת-גרף זה:

$h_{i}:={\begin{cases}1/{\sqrt {|A_{j}|}}&{\mbox{if}}\ v_{i}{\mbox{ in }}A_{j}\\0&{\mbox{otherwise}}\end{cases}}$ מתקיים השוויון הבא: $h'Lh={\frac {cut(A_{j},{\overline {A_{j}}})}{|A_{j}|}}$

ובאופן כללי, עבור מטריצה $H$ מגודל $n\times k$ אשר עמודותיה הן וקטורים אינדיקטורים נקבל את השוויון הבא: $Trace(H'LH)=\sum _{i\,=\,1}^{k}{\frac {cut(A_{i}{\overline {A_{i}}})}{|A_{i}|}}$ כלומר, מציאת מינימום לפונקציה $RatioCut$ שקולה למציאת מינימום לביטוי $Trace(H'LH)$ תחת ההנחה שעמודות $H$ הן וקטורים אינדיקטורים וכן שהן אורתונורמליות זו לזו.

בעיה זו אמנם עדיין NP-קשה, אך ניתן למצוא פתרון מקורב לבעיה זו אם מסירים את הדרישה שהעמודות הן אינדיקטורים, ומניחים רק שהן אורתונורמליות, וזוהי כבר בעיה ידועה אשר פתרונה נתון על ידי הווקטורים העצמיים המתאימים ל- $k$ הערכים העצמיים הקטנים ביותר של המטריצה $L$ , וסיבוכיות החישוב שלה פולינמיאלית. באופן דומה ניתן למצוא את הווקטורים העצמיים המתאימים ל- $k$ הערכים העצמיים הקטנים ביותר של מטריצת הלפלסיאן המנורמלת $L^{rw}$ כדי לקרב את פונקציית המטרה $nCut$ ^[5].

ניתן להבין קירוב זה בצורה אינטואיטיבית - קל להראות כי עבור גרף בעל $k$ רכיבי קשירות, הפתרון האופטימלי עבור $k$ אשכולות יהיה חלוקת הגרף לאותם $k$ רכיבי קשירות, והווקטורים העצמיים המתאימים לערך העצמי 0 הם אינדיקטורים של רכיבי הקשירות^[6], כמו כן אלו הם הערכים העצמיים הקטנים ביותר, מאחר שידוע כי הערכים העצמיים של מטריצת הלפלסיאן אינם שליליים. כעת עבור גרף דומה המכיל בנוסף "רעש קטן", כלומר קשתות מעטות בין רכיב קשירות אחד למשנהו, נוכל לצפות שהערכים העצמיים ישתנו מעט, והמרחב הווקטורי הנפרש על ידי הווקטורים העצמיים המתאימים להם ישתנה מעט אף הוא, כלומר ניתן להבין כי בהינתן רעש קטן, קירוב הפתרון יהיה בסביבת הפתרון המקורי^[7].

חשוב לציין כי אין ערובה לכך שקירוב זה אכן יניב את הפתרון לבעיה המקורית (המניחה שהווקטורים הם אינדיקטורים), דוגמאות קלאסיות לכך נסקרו בספרות במהלך השנים^[8]. ובנוסף, בעוד פתרון הבעיה המקורית הוא אוסף אינדיקטורים המחלק את המידע ל $k$ קבוצות בצורה מפורשת, קירוב הפתרון כבר מכריח אותנו לנתח את הווקטורים העצמיים בעזרת אלגוריתם נוסף (למשל בעזרת אלגוריתם k-מרכזים).

דוגמאות

להלן דוגמה לגרף בעלי 2 רכיבי קשירות, וכן גרף זהה עם תוספת של רעש ההופך אותו לבעל רכיב קשירות אחד,

גרף	מטריצת לפלסיאן	2 וקטורים עצמיים ראשונים
	$\left({\begin{array}{rrrrrr}20&-20&0&0\\-20&20&0&0\\0&0&30&-30\\0&0&-30&30\\\end{array}}\right)$	$\left({\begin{array}{rrrr}0\\0\\1\\1\end{array}}\right)$ $\left({\begin{array}{rrrr}1\\1\\0\\0\end{array}}\right)$

גרף	מטריצת לפלסיאן	2 וקטורים עצמיים ראשונים
	$\left({\begin{array}{rrrrrr}21&-20&0&-1\\-20&20&0&0\\0&0&30&-30\\-1&0&-30&31\\\end{array}}\right)$	$\left({\begin{array}{rrrr}-0.991\\-1.042\\1.034\\1\end{array}}\right)$ $\left({\begin{array}{rrrr}1\\1\\1\\1\end{array}}\right)$

ניתן לראות בגרף השני את הקשת המחברת את הצומת הראשון והשלישי כרעש קטן, ואכן הווקטורים העצמיים פורשים מרחב קרוב למרחב הנפרש על ידי האינדיקטורים של "רכיבי הקשירות" שכן מתקיים: ${\frac {v_{1}+v_{2}}{2}}=\left({\begin{array}{rrrr}0.0045\\-0.021\\1.017\\1\end{array}}\right),{\frac {v_{1}-v_{2}}{2}}=\left({\begin{array}{rrrr}0.9955\\1.021\\-0.017\\0\end{array}}\right)$

הווקטור העצמי הראשון

חשוב לציין כי במרבית המקרים הקלט המתקבל הוא גרף קשיר, ומכאן 0 הוא ערך עצמי מריבוי 1, והווקטור העצמי המתאים לו יהיה וקטור אינדיקטור של כל צמתי הגרף, כלומר - $v=(1,1,...,1)$ .

שימוש בוקטור זה לא מוסיף מידע על הקשר בין צמתי הגרף, ולכן ניתן לא להשתמש בו כאשר מחלקים את הגרף לאשכולות. לדוגמה, עבור חלוקה לשני אשכולות, ניתן להסתכל ישירות על הווקטור העצמי השני הקטן ביותר (הנקרא Fiedler Vector), ולחלק את הגרף על פיו, אשכול אחד יהיה קבוצת הצמתים בעלי ערך חיובי בוקטור, והאשכול השני יהיו קבוצת הצמתים בעלי ערך שלילי.

אלגוריתם בסיסי

כפי שהוזכר, בגרף בעל $k$ רכיבי קשירות, הפתרון המתקבל עבור $k$ הווקטורים העצמיים הראשונים פורש את אותו המרחב הנפרש גם כן על ידי האינדיקטורים של רכיבי הקשירות, ועל כן מציאת $k$ האשכולות תתקבל על ידי יצירת מטריצה בגודל $k\times n$ מהווקטורים העצמיים, והבאתה לצורתה הקנונית.

אמנם בדרך כלל הגרף יכיל רעש, או שמראש מספר האשכולות לא יהיה ידוע, ועל כן מספר רכיבי הקשירות יהיה שונה ממספר האשכולות ויידרש עיבוד נוסף של המידע, דרך נפוצה אחת היא שימוש באלגוריתם k-מרכזים. בנוסף, במקרים רבים קלט המערכת אינו מטריצת לפלסיאן אלא, למשל, אוסף נקודות במימד גבוה כלשהו, במקרה כזה נדרש להגדיר פונקציית סמיכות המציינת את עוצמת הקשר בין כל שני איברים כדי ליצור את מטריצת הלפלסיאן. האלגוריתם השלם יכלול את השלבים הבאים:

יצירת גרף ומטריצת לפלסיאן מהמידע
חישוב $k$ הווקטורים העצמיים הראשונים של המטריצה
הטלת מטריצת הווקטורים הללו למרחב $R^{k}$ , השורה ה $l$ של מטריצה זו היא הטלת הצומת $l$ בגרף אל $R^{k}$
הפעלת אלגוריתם לניתוח אשכולות במרחב זה (לדוגמה אלגוריתם k-מרכזים).

להלן קטע קוד המממש את האלגוריתם המתואר^[9]

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_circles, make_blobs
from sklearn.cluster import KMeans
from sklearn.neighbors import NearestNeighbors
from scipy.linalg import eigh

# Step 0: Generate a dataset (Circles + Blobs)
np.random.seed(42)
X_circles, y_circles = make_circles(n_samples=400, factor=0.5, noise=0.05)
X_blobs, y_blobs = make_blobs(n_samples=100, centers=[(-1.5, -1.5)], cluster_std=0.1)
X = np.vstack((X_circles, X_blobs))
y_true = np.hstack((y_circles, y_blobs + 2))  # Adjust labels for clarity

# Step 1: Create the graph Laplacian matrix using epsilon-neighborhood
epsilon = 0.2
nbrs = NearestNeighbors(radius=epsilon).fit(X)
adj_matrix = nbrs.radius_neighbors_graph(X, mode='connectivity').toarray()
D = np.diag(np.sum(adj_matrix, axis=1))
L = D - adj_matrix

# Step 2: Compute the eigenvectors of the Laplacian matrix
eigenvalues, eigenvectors = eigh(L)
eigenvectors = eigenvectors[:, 0:3]

# Step 3+4: Analyze the data with k-means
kmeans = KMeans(n_clusters=3, random_state=42)
labels = kmeans.fit_predict(eigenvectors)

# Draw the original labels and clustering labels
plt.figure(figsize=(14, 6))

plt.subplot(1, 2, 1)
plt.scatter(X[:, 0], X[:, 1], c=y_true)
plt.title('Original Data (True Labels)')


plt.subplot(1, 2, 2)
plt.scatter(X[:, 0], X[:, 1], c=labels)
plt.title('Spectral Clustering Result')


plt.tight_layout()
plt.show()

מעבר מדאטה במימד גבוה למטריצת שכנות

קיימות מספר שיטות ליצור גרף בהינתן אוסף נקודות, הבולטות שבהן הן:

גרסה של אלגוריתם שכן קרוב - נגדיר קשת $U_{ij}$ $U_{ij}$ אם $i$ $i$ הוא אחד מ $k$ $k$ השכנים הקרובים של $j$ $j$ , חשוב לציין כי:
- הערך $k$ נתון לבחירה, ככל שיהיה גדול יותר כך הגרף יכיל יותר קשתות.
- הגרף המתקבל בדרך זו הוא מכוון, וזו תופעה שאיננה רצויה, לכן ניתן לבחור האם להגדיר קשת רק במקרים שבהם שתי הקשתות $U_{ij}$ , $U_{ji}$ קיימות (אלגוריתם שכן קרוב הדדי), או גם במקרים שבהם לפחות אחת משתי הקשתות קיימות.
- הגרף המתקבל איננו ממושקל
סביבת $\epsilon$ $\epsilon$ - נגדיר קשת $U_{ij}$ $U_{ij}$ אם המרחק בין שני הצמתים קטן או שווה ל- $\epsilon$ $\epsilon$ , נשים לב כי:
- הערך $\epsilon$ נתון לבחירה, ככל שיהיה גדול יותר כך הגרף יכיל יותר קשתות.
- הגרף המתקבל בדרך זו איננו מכוון ואיננו ממושקל.
שימוש במטריקת מרחק המודדת קרבה בין שתי נקודות, פונקציה שכיחה בהקשר זה היא - $s(x_{i},x_{j})=e^{\frac {-||x_{i}-x_{j}||^{2}}{2\sigma ^{2}}}$ $s(x_{i},x_{j})=e^{\frac {-||x_{i}-x_{j}||^{2}}{2\sigma ^{2}}}$ , בשימוש במטריקה זו:
- הערך $\sigma$ נתון לבחירה, ככל שיהיה גדול יותר כך נקודות בסביבה רחוקה יותר יחוברו בקשת בעלת משקל גבוה.
- הגרף המתקבל הוא ממושקל ואיננו מכוון.
- הגרף המתקבל הוא שלם, ועל כן מטריצת הלפלסיאן לא תהיה דלילה, וזמני החישוב יהיו ארוכים יותר.

מאפיינים מרכזיים בהשוואה לשיטות שונות

סוג הקלט

חשוב לשים לב כי האלגוריתם מתבסס על דמיון וקירבה בין איברים, ולכן ניתן להשתמש בקלט בעל אופי מגוון. בניגוד לאלגוריתם k-מרכזים או אלגוריתם DBSCAN(אנ') לדוגמה, אשר מתאימים יותר לקלט ממימד $R^{n}$ , כאן ניתן להשתמש גם בערכים בוליאנים (כאשר ניתן להגדיר מטריקת מרחק המסתמכת על מרחק ג'אקרד לדוגמה).

אופן החישוב

השלבים ה"כבדים" הם יצירת מטריצת הלפלסיאן, וחישוב הווקטורים העצמיים. עבור $n$ נקודות קלט, בניית מטריצת לפלסיאן מלאה פרופורציונלית למספר הקשתות - כלומר סיבוכיות של ${\mathcal {O}}(n^{2})$ , ומספר פעולות החישוב עבור כל הווקטורים העצמיים במדויק פרופורציונלי ל $n^{3}$ .

חשוב לציין כי לשימוש מעשי נדרש רק קירוב של הווקטורים העצמיים, ובנוסף ביישומים רבים, תחת הנחות ריאליות של ידע מקדים על הקלט ניתן להגביל את מספר הקשתות של כל צומת, ובכך עבור $k$ וקטורים עצמיים ניתן להוריד את הסיבוכיות עד לסדר גודל של ${\mathcal {O}}(k*n^{1.5})$ בעזרת שימוש באלגוריתם לנצוש^[10]^[11]. לכן זמן הריצה תלוי מאוד באופי הקלט ואופן מימוש האלגוריתם.

עוד יש לציין כי מספר הממדים של הקלט איננו משפיע על זמני הריצה של שלבים אלה בצורה ישירה, מה שמבדיל אותו מאלגוריתמים רבים אשר יעילותם פוחתת וזמני הריצה שלהם עולים משמעותית תחת קלט של מידע ממימד גבוה.

התפלגות הקלט

אחד היתרונות המרכזיים של האלגוריתם הוא העובדה שאין הנחות יסוד על האשכולות. בניגוד ל- $k$ -מרכזים המניח שכל אשכול הוא בעל צורה כללית של ספירה, או DBSCAN המניח צפיפות אחידה באשכול, כאן פונקציית המטרה לא מתייחסת בצורה ישירה לתכונות גלובליות אלה של האשכול.

ראו גם

קישורים חיצוניים

מאמר מקיף מאת פרופסור ulrike von luxburg.
סדרת סרטונים בנושא מאוניברסיטת סטנפורד
תרגול קוד פייתון בסיסי עבור סגמנטציה מאת scikit-learn
תיעוד ספריית sklearn.cluster.SpectralClustering

הערות שוליים

^ Wagner, D. and Wagner, F. (1993). Between min cut and graph bisection
^ Hagen and Kahng, 1992
^ Shi and Malik, 2000
^ להוכחה המלאה ראו A Tutorial on Spectral Clustering by Ulrike von Luxburg בפרק 5.2
^ להוכחה המלאה ראו A Tutorial on Spectral Clustering by Ulrike von Luxburg בפרק 5.3
^ ליתר דיוק יש לומר כי המרחב הנפרש על ידי הווקטורים העצמיים זהה למרחב הנפרש על ידי אינדיקטורים, מאחר שלא ניתן לשייך וקטור עצמי אחד מסוים אל ערך עצמי בעל ריבוי גאומטרי גדול מאחד אלא יש לדבר על המרחב הווקטורי המתאים לערך העצמי
^ davis and kahan, 1970
^ ראו לדוגמה Guattery and Miller, 1978
^ קוד זה עוקב אחר השלבים המרכזיים של האלגוריתם התאורטי, ברמה המעשית ניתן להשתמש בספריות ייעודיות כגון sklearn.cluster.SpectralClustering
^ למעשה השילוב של הסתפקות בקירוב הווקטורים העצמיים עם הגבלת מספר הקשתות המחוברות לכל צומת הוא שהופך את השימוש באלגוריתם לנצוש לאטרקטיבי מאחר שזמן החישוב של אלגוריתם זה יורד משמעותית עבור מטריצות דלילות, בהיעדר הנחה זו ניתן להגיע לסיבוכיות של $k*n^{2}$ עבור $k$ הווקטורים העצמיים הראשונים
^ ראו בהרחבה בshi and malik, 2000, פרק 3.1

[1] Wagner, D. and Wagner, F. (1993). Between min cut and graph bisection

[2] Hagen and Kahng, 1992

[3] Shi and Malik, 2000

[4] להוכחה המלאה ראו A Tutorial on Spectral Clustering by Ulrike von Luxburg בפרק 5.2

[5] להוכחה המלאה ראו A Tutorial on Spectral Clustering by Ulrike von Luxburg בפרק 5.3

[6] ליתר דיוק יש לומר כי המרחב הנפרש על ידי הווקטורים העצמיים זהה למרחב הנפרש על ידי אינדיקטורים, מאחר שלא ניתן לשייך וקטור עצמי אחד מסוים אל ערך עצמי בעל ריבוי גאומטרי גדול מאחד אלא יש לדבר על המרחב הווקטורי המתאים לערך העצמי

[7] vis and kahan, 1970

[8] ראו לדוגמה Guattery and Miller, 1978

[9] קוד זה עוקב אחר השלבים המרכזיים של האלגוריתם התאורטי, ברמה המעשית ניתן להשתמש בספריות ייעודיות כגון sklearn.cluster.SpectralClustering

[10] למעשה השילוב של הסתפקות בקירוב הווקטורים העצמיים עם הגבלת מספר הקשתות המחוברות לכל צומת הוא שהופך את השימוש באלגוריתם לנצוש לאטרקטיבי מאחר שזמן החישוב של אלגוריתם זה יורד משמעותית עבור מטריצות דלילות, בהיעדר הנחה זו ניתן להגיע לסיבוכיות של $k*n^{2}$ עבור $k$ הווקטורים העצמיים הראשונים

[11] ראו בהרחבה בshi and malik, 2000, פרק 3.1

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]