ג.7 רגרסיה לינארית — משקל אחד

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

זו בדיוק המטרה של רגרסיה לינארית — Linear Regression: למצוא קו ישר שמתאר היטב את הקשר בין הקלט x לבין התוצאה y. הקו מתאר את מגמת הנתונים, ואינו חייב לעבור בכל הנקודות. ברגע שיש בידינו הקו, אפשר להציב בו כל ערך של x — גם ערך שלא הופיע בנתונים — ולקבל הערכה של y. היכולת הזו, להסיק מהדוגמאות שראינו על מקרים שלא ראינו, נקראת הכללה — Generalization, והיא לב־לבה של למידת מכונה: המודל אינו לומד את הנתונים בעל־פה, אלא לומד את הכלל שמסתתר מאחוריהם.

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

בפרק זה נתחיל במקרה הפשוט ביותר: ישר העובר בראשית הצירים, שהמשוואה שלו היא y = wx; המשקל w הוא השיפוע, והוא הפרמטר היחיד שנלמד. תחילה נתאים את הישר לנקודה אחת, ואחר כך לכמה נקודות יחד. בפרק הבא נוסיף לישר פרמטר שני, ההטיה, שיאפשר לו שלא לעבור בראשית.

השיעור וההרצאות באתר של גלעד מרקמן

חומרי הליווי: 5. רגרסיה לינארית (עותק מקומי) · מחברת רגרסיה לינארית (עותק מקומי)

נייבא את הספריות שישמשו אותנו לאורך הפרק: PyTorch לחישובים ולאימון, NumPy ליצירת טווחי ערכים, ו־Matplotlib לציור הנקודות והישרים:

import torch
import numpy as np
import matplotlib.pyplot as plt

רגרסיה בנקודה אחת

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

הכנת הנתונים

לפני שנאמן דבר, נראה מה קורה עם ניחוש שרירותי. נציג תחילה ישר בעל שיפוע 1. הוא מנבא עבור x=2 את הערך 2, בעוד שהתוצאה הרצויה היא 4.5. הפער בין השניים הוא הטעות שנרצה להקטין:

x, y = 2, 4.5
w = 1

y_pre = w * x

plt.scatter([0,x],[0,y])
plt.plot([0,x], [0,y_pre], color = 'red')
plt.title("y_predict = w * x")
plt.ylabel("y_predict")
plt.xlabel("x")

פונקציית ההפסד

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

$$ loss=(y_{predict}-y)^2=(wx-y)^2 $$

הריבוע דואג שההפסד יהיה תמיד חיובי, בין שהתחזית גבוהה מדי ובין שהיא נמוכה מדי, ושטעויות גדולות "ייענשו" יותר מטעויות קטנות. שימו לב לשינוי נקודת המבט: בגרף הקודם ציירנו את y כתלות ב־x, ועכשיו x ו־y קבועים ואנו שואלים כיצד ההפסד משתנה כשמשנים את w. נצייר את ההפסד עבור שיפועים בין ‎−4 ל־9. כעת הציר האופקי הוא המשקל w, והציר האנכי הוא ההפסד.

w_np = np.linspace(-4, 9, 100)

# loss_np = (w_np*2 - 4.5)**2
loss_np = (w_np*x - y)**2
plt.plot(w_np, loss_np)
plt.title("loss = (w * x - y)^2;   x = 2; y = 4.5")
plt.ylabel("loss")
plt.xlabel("W")

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

הגדרת המודל

מודל — Model הוא הפונקציה שמקבלת קלט ומחזירה תחזית; בהמשך הספר המודלים יהיו רשתות נוירונים, אך כאן המודל שלנו פשוט מאוד: הוא מחשב תחזית באמצעות הכפלת הקלט במשקל שנלמד. נגדיר את הנתונים כטנסורים, כדי ש־Autograd יוכל לחשב נגזרות. את המשקל שנאמן נאתחל ל־‎−1.5 עם requires_grad=True — ערך התחלתי רחוק בכוונה מהפתרון, כדי שנוכל לראות את האימון מתקדם — ונבחר קצב למידה 0.01:

X = torch.tensor(2)
Y = torch.tensor(4.5)
W =  torch.tensor(-1.5, dtype=torch.float32, requires_grad=True)
learning_rate = 0.01

def model (X):
    return W * X

פונקציית ההפסד והאופטימייזר

נתרגם עתה לקוד את שני המרכיבים הנותרים. הפונקציה loss מחשבת את ריבוע השגיאה שהגדרנו לעיל. האופטימייזר, כפי שראינו בפרק ג.5, הוא הרכיב שמבצע את צעד העדכון של Gradient Descent; הוא מקבל את רשימת הפרמטרים שעליו לעדכן — כאן רק W — ואת קצב הלמידה, ובכל צעד מזיז את W בכיוון ההפוך לנגזרת:

def loss (Y_predict, Y):
    return (Y_predict - Y)**2  #(W * X - Y)**2

optim = torch.optim.SGD([W], lr=learning_rate) # W = W - grad * learning_rate

לולאת האימון

לולאת האימון היא אותה לולאה שהכרנו בפרקי Gradient Descent, אלא שהפעם הפונקציה שאנו ממזערים היא ההפסד של המודל. כל מעבר על הלולאה נקרא אפוק — epoch. בכל איטרציה מחשבים תחזית והפסד (המעבר קדימה), מחשבים את הנגזרת (המעבר לאחור), מעדכנים את המשקל ומאפסים את הנגזרת. אחת לעשר איטרציות נדפיס את המצב כדי לעקוב אחר ההתקדמות:

epochs = 300
for epoch in range(epochs):
    # Forward
    Y_predict = model(X)

    # Backward - calculate new grads
    l = loss(Y_predict, Y)
    l.backward()

    if (epoch) % 10 == 0:
        print (f"epoch={epoch} W={W:.3f}, Y_predict={Y_predict:.3f}, loss={l:.6f} grad={W.grad}")

    # Update weight
    optim.step()  # w = w - grad * LR

    # zero grads
    optim.zero_grad()

שורת הפלט הראשונה והאחרונה של הלולאה:

פלט

epoch=0 W=-1.500, Y_predict=-3.000, loss=56.250000 grad=-30.0
epoch=290 W=2.250, Y_predict=4.500, loss=0.000000 grad=-1.1444091796875e-05

נקרא את הפלט. בתחילת האימון המשקל הוא ‎−1.5, התחזית היא ‎−3 (רחוקה מאוד מ־4.5), ההפסד 56.25 והנגזרת ‎−30: נגזרת שלילית גדולה פירושה שיש להגדיל את המשקל בצעדים גדולים. בסוף האימון המשקל הגיע ל־2.25, התחזית היא בדיוק 4.5, ההפסד אפס והנגזרת כמעט אפס — הגענו לתחתית הפרבולה שראינו בגרף ההפסד.

הצגת התוצאה

נחשב שוב את התחזית ואת ההפסד ונצייר את הישר שהתקבל:

Y_predict = model(X)
l = loss(Y_predict, Y)
print (f"W={W:.3f}, Y_predict={Y_predict:.3f}, loss={l:.6f}")

plt.scatter(x,y)
plt.plot([0,x], [0,Y_predict.item()], color = 'red')
# plt.plot([0, 2], [0, 2.4], color='b')
plt.show()

השיפוע שהתקבל הוא בקירוב 2.25, ולכן התחזית עבור x=2 היא בקירוב 4.5. הישר "למד" את הנקודה בלי שאמרנו לו את התשובה — מצאנו אותה דרך הקטנת ההפסד בלבד.

רגרסיה עם כמה נקודות

נקודה אחת היא מקרה מלאכותי: תמיד אפשר להעביר ישר דרך הראשית ודרך נקודה אחת, וההפסד יורד לאפס. המצב המעניין באמת, זה שתיארנו בפתיחת הפרק, הוא כשיש נקודות רבות שאינן על קו אחד, ואנו מחפשים ישר יחיד שמתאר את כולן "הכי טוב שאפשר". כעת נחפש משקל אחד שמתאים לכמה נקודות יחד. נכין ארבע נקודות. ערכי y מתקבלים מהכפלת x ב־2 ומהוספת מספר אקראי בין 0 ל־1 — הרעש האקראי מדמה את הגורמים שאינם בשליטתנו, כמו בדוגמת המכירות, ובגללו שום ישר לא יעבור בכל הנקודות. הפקודה torch.manual_seed(1) קובעת את המספרים האקראיים כך שנקבל אותם ערכים בכל הרצה:

torch.manual_seed(1)
X = torch.tensor([0, 1, 2, 3], dtype=torch.float32)
Y = X * 2 + torch.rand(4,)
print (Y)
plt.plot (0,0)
plt.scatter(X,Y)
plt.plot([0,4], [0,4])

פלט

tensor([0.7576, 2.2793, 4.4031, 6.7347])

כשיש כמה נקודות, לכל אחת מהן טעות משלה, ואנו זקוקים למספר אחד שמסכם את כולן. הדרך המקובלת היא לחשב את ממוצע ריבועי השגיאות — MSE, Mean Squared Error:

$$ loss=\frac{1}{n}\sum_{i=1}^{n}(wx_i-y_i)^2 $$

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

גרף פונקציית ההפסד

כמו במקרה של נקודה אחת, נצייר תחילה את פונקציית ההפסד כדי לראות מה אנו עומדים למזער. גרף ההפסד מאפשר לראות כיצד איכות ההתאמה משתנה עם המשקל. נחשב את ההפסד עבור 100 משקלים בין ‎−4 ל־9. במקום לולאה על המשקלים נשתמש בכפל מטריצות: W@X יוצר טבלת תחזיות: שורה לכל משקל ועמודה לכל נקודה; mean(1) מחשבת את הממוצע לאורך כל שורה, כלומר את ההפסד הממוצע לכל משקל. השורות המוערות בקוד מראות את החישוב המקביל בלולאה רגילה.

W = torch.tensor(np.linspace(-4, 9, 100), dtype=torch.float32)
W = W.reshape(-1,1)
X = X.reshape(1,-1)
print(W.shape, X.shape)
loss = ((W@X-Y)**2).mean(1)
print(loss)
# print ("W", W)
# print ("X", X)
# print ("W@X",W@X)
# loss = torch.zeros(100)
# for i in range(W.size()[0]):
#   loss[i] = ((W[i]*X - Y)**2).mean()

plt.plot(W.numpy(), loss.numpy())
plt.title("Loss = ((w * x - Y)**2).mean()")
plt.ylabel("Loss")
plt.xlabel("W")

פלט

torch.Size([100, 1]) torch.Size([1, 4])
tensor([1.3621e+02, 1.3054e+02, 1.2499e+02, 1.1956e+02, 1.1425e+02, 1.0906e+02,
        1.0399e+02, 9.9046e+01, 9.4220e+01, 8.9515e+01, 8.4931e+01, 8.0467e+01,
        7.6124e+01, 7.1902e+01, 6.7800e+01, 6.3819e+01, 5.9959e+01, 5.6219e+01,
        5.2600e+01, 4.9102e+01, 4.5725e+01, 4.2468e+01, 3.9332e+01, 3.6317e+01,
        3.3422e+01, 3.0648e+01, 2.7995e+01, 2.5462e+01, 2.3051e+01, 2.0760e+01,
        1.8589e+01, 1.6539e+01, 1.4610e+01, 1.2802e+01, 1.1114e+01, 9.5475e+00,
        8.1013e+00, 6.7758e+00, 5.5710e+00, 4.4869e+00, 3.5235e+00, 2.6808e+00,
        1.9588e+00, 1.3575e+00, 8.7690e-01, 5.1701e-01, 2.7782e-01, 1.5933e-01,
        1.6155e-01, 2.8446e-01, 5.2808e-01, 8.9240e-01, 1.3774e+00, 1.9832e+00,
        2.7096e+00, 3.5567e+00, 4.5245e+00, 5.6131e+00, 6.8223e+00, 8.1522e+00,
        9.6029e+00, 1.1174e+01, 1.2866e+01, 1.4679e+01, 1.6612e+01, 1.8667e+01,
        2.0841e+01, 2.3137e+01, 2.5553e+01, 2.8090e+01, 3.0748e+01, 3.3526e+01,
        3.6425e+01, 3.9445e+01, 4.2586e+01, 4.5847e+01, 4.9229e+01, 5.2731e+01,
        5.6354e+01, 6.0098e+01, 6.3963e+01, 6.7948e+01, 7.2055e+01, 7.6281e+01,
        8.0629e+01, 8.5097e+01, 8.9686e+01, 9.4395e+01, 9.9226e+01, 1.0418e+02,
        1.0925e+02, 1.1444e+02, 1.1975e+02, 1.2519e+02, 1.3074e+02, 1.3642e+02,
        1.4221e+02, 1.4813e+02, 1.5417e+02, 1.6033e+02])

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

הכנת הנתונים והמודל לאימון

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

torch.manual_seed(1)
X = torch.tensor([0, 1, 2, 3, 4], dtype=torch.float32)
Y = X * 2 + torch.rand(5,)
W = torch.tensor(8.0, requires_grad=True)
learning_rate = 0.01
print (X, Y, W)

def Model (X):
    return W * X
Y_predict = Model(X)
print (Y_predict)

פלט

tensor([0., 1., 2., 3., 4.]) tensor([0.7576, 2.2793, 4.4031, 6.7347, 8.0293]) tensor(8., requires_grad=True)
tensor([ 0.,  8., 16., 24., 32.], grad_fn=<MulBackward0>)

התחזיות ההתחלתיות (0, 8, 16, 24, 32) גדולות בהרבה מהערכים האמיתיים (כ־0.76 עד 8.03), כצפוי ממשקל 8 כשהשיפוע האמיתי הוא בערך 2.

הגדרת ההפסד והאופטימייזר

ההבדל היחיד לעומת המקרה של נקודה אחת הוא בפונקציית ההפסד: במקום ריבוע שגיאה יחיד נחשב את ה־MSE, כלומר נוסיף .mean() על טנסור השגיאות. נגדיר את ממוצע ריבועי השגיאות כפונקציה וניצור אופטימייזר עבור המשקל:

def Loss (Y_predict, Y):
    return ((Y_predict - Y)**2).mean() #((W * X - Y)**2).mean()

optimizer = torch.optim.SGD([W], lr=learning_rate) # W = W - W.grad * LR

אימון המודל

לולאת האימון זהה לחלוטין לזו של הנקודה היחידה — זהו יתרון גדול של הגישה: אותו קוד עובד לנקודה אחת, לחמש נקודות ולאלפי נקודות. נבצע 200 איטרציות על חמש הנקודות:

for epoch in range(200):

    # forward
    Y_predict = Model(X)

    # backward
    loss = Loss(Y_predict, Y)
    loss.backward()

    if epoch % 1 == 0:
        print(f"epoch= {epoch} W= {W.item():.3f} loss={loss.item():.3f} grad= {W.grad.item():.3f}")

    # Update weight
    optimizer.step()

    # zero grads
    optimizer.zero_grad()

שורת הפלט הראשונה והאחרונה:

פלט

epoch= 0 W= 8.000 loss=208.095 grad= 70.637
epoch= 199 W= 2.114 loss=0.194 grad= 0.000

בתחילת האימון ההפסד הוא 208 והנגזרת חיובית וגדולה (כ־70.6), ולכן המשקל יקטן במהירות. אחרי 200 אפוקים הנגזרת התאפסה, כלומר הגענו למינימום, והמשקל התייצב על 2.114 — קרוב ל־2, השיפוע שממנו יצרנו את הנתונים.

הצגת ההתאמה

נציג את המשקל ואת הישר שהתקבל. הבלוק with torch.no_grad() אומר ל־PyTorch שלא לעקוב אחר החישוב לצורך נגזרות, כי כאן אנו רק מציירים ולא מאמנים:

print(W)
print(f"END: W= {W.item():.3f} loss={loss.item():.3f}")
plt.scatter(X, Y, color='b')
plt.xlabel("x")
plt.ylabel("Y")
with torch.no_grad():
    plt.plot(X, Model(X), color='r')
plt.show()

פלט

tensor(2.1136, requires_grad=True)
END: W= 2.114 loss=0.194

ההפסד אינו אפס, משום שהישר אינו עובר בכל הנקודות בדיוק. זה אינו כישלון אלא בדיוק מה שרצינו: הישר תופס את המגמה הכללית ולא את הרעש האקראי של כל נקודה. עכשיו אפשר לעשות את מה שדיברנו עליו בפתיחה — להציב בישר ערך x חדש, למשל 2.5, ולקבל הערכה סבירה ל־y, אף שמעולם לא ראינו נקודה כזו. זוהי ההכללה.

לסיום, נכין את הנתונים לפרק הבא. הפעם ערכי y מתקבלים מהנוסחה y = 2x + 5, ולכן הנקודות מונחות על ישר שאינו עובר בראשית הצירים. המודל שלנו, y = wx, אינו מסוגל לתאר ישר כזה ולו רק מפני שהוא נאלץ לעבור דרך (0, 0); כדי לפתור זאת נוסיף בפרק הבא פרמטר שני — ההטיה:

x = np.array([0, 1, 2, 3, 4], dtype=float)
y = x * 2 + 5

print(y)

פלט

[ 5.  7.  9. 11. 13.]