Gauss và tiểu hành tinh biến mất: phương pháp bình phương nhỏ nhất ra đời thế nào

Một ngôi sao không đứng yên

Đêm đầu tiên của thế kỷ 19, ngày 1 tháng 1 năm 1801, nhà thiên văn Giuseppe Piazzi ở đài quan sát Palermo, đảo Sicilia, đang lập bản đồ các ngôi sao. Ông ghi lại vị trí của một đốm sáng mờ trong chòm Kim Ngưu. Đêm hôm sau, đốm sáng ấy đã dịch đi một chút. Đêm hôm sau nữa, nó lại dịch tiếp.

Sao thì không di chuyển như thế. Piazzi đã tìm thấy một thiên thể mới trong Hệ Mặt Trời. Ông đặt tên nó là Ceres, theo tên nữ thần nông nghiệp La Mã, vị thần bảo hộ của Sicilia.

Piazzi theo dõi Ceres được 24 lần trong khoảng 41 ngày. Ngày 11 tháng 2, ông lâm bệnh và phải ngừng quan sát. Khi ông khỏe lại, Ceres đã trôi vào vùng trời gần Mặt Trời, bị chìm trong ánh sáng ban ngày[2]. Nó sẽ không thể được nhìn thấy trong nhiều tháng.

Vấn đề là trong 41 ngày ấy, Ceres chỉ vạch một cung rất ngắn trên bầu trời, chưa tới mười độ. Từ một mẩu đường cong nhỏ như vậy, phải đoán ra cả một quỹ đạo quanh Mặt Trời, rồi tính xem nhiều tháng sau nó sẽ hiện ra ở đâu giữa bầu trời đầy sao. Không ai biết cách làm việc đó. Các nhà thiên văn châu Âu đưa ra những dự đoán khác nhau, và bầu trời ở những chỗ họ chỉ vẫn trống trơn.

Thiên thể mới được phát hiện ấy có nguy cơ bị lạc mất mãi mãi.

Cuối năm 1801, một chàng trai 24 tuổi ở Brunswick, nước Đức, chưa từng làm thiên văn chuyên nghiệp, gửi đi một dự đoán khác hẳn tất cả những người còn lại. Và Ceres hiện ra gần đúng ở chỗ anh chỉ.

Đây là tập thứ ba trong loạt bài về các nhà toán học. Cuối bài, bạn sẽ thấy ý tưởng giúp tìm lại Ceres là gì, và vì sao nó đang chạy trong chiếc điện thoại của bạn mỗi khi bạn mở bản đồ.

Cảnh sát bầu trời

Phát hiện của Piazzi không đến từ hư không. Từ cuối thế kỷ 18, nhiều nhà thiên văn tin rằng giữa quỹ đạo của sao Hỏa và sao Mộc phải có một hành tinh. Khoảng cách giữa hai quỹ đạo ấy rộng bất thường, và một quy luật số học khi đó đang rất được tin tưởng, quy tắc Titius–Bode, “dự đoán” có một hành tinh ở đúng vị trí ấy.

Năm 1800, nhà thiên văn Franz Xaver von Zach tập hợp một nhóm đồng nghiệp khắp châu Âu, chia bầu trời thành từng vùng để mỗi người dò tìm một phần. Họ tự gọi mình là “cảnh sát bầu trời”. Trớ trêu thay, Piazzi chưa kịp nhận được lời mời tham gia thì đã tìm ra Ceres một cách tình cờ.

Ceres nằm đúng ở khoảng cách mà mọi người chờ đợi. Ngày nay ta biết nó không phải là một hành tinh lớn, mà là thiên thể lớn nhất trong vành đai tiểu hành tinh, đường kính khoảng 940 km, và từ năm 2006 được xếp vào nhóm hành tinh lùn. Nhưng vào năm 1801, nó là “hành tinh bị thất lạc”, và cả giới thiên văn châu Âu muốn tìm lại nó.

Chàng trai ở Brunswick

Carl Friedrich Gauss sinh năm 1777 trong một gia đình nghèo ở Brunswick. Cha ông là thợ xây kiêm thợ làm vườn. Tài năng toán học của Gauss lộ ra rất sớm, sớm tới mức công tước xứ Brunswick đã cấp học bổng cho ông học đại học.

Năm 1796, khi chưa tròn 19 tuổi, Gauss chứng minh rằng có thể dựng một đa giác đều 17 cạnh chỉ bằng thước kẻ và compa. Đây là câu hỏi đã bị bỏ ngỏ từ thời Hy Lạp cổ đại, suốt hơn hai nghìn năm. Ông vui đến mức quyết định theo toán học thay vì ngôn ngữ cổ điển, và về sau muốn hình đa giác 17 cạnh được khắc trên bia mộ mình. Năm 1801, cũng là năm Ceres biến mất, ông xuất bản cuốn Disquisitiones Arithmeticae, cuốn sách đặt nền móng cho lý thuyết số hiện đại[1].

Tức là vào năm 1801, Gauss là một nhà lý thuyết số trẻ, không phải nhà thiên văn. Nhưng bài toán Ceres hấp dẫn ông vì một lý do: đó là bài toán tính toán, và người ta đang thiếu một phương pháp. Gauss bắt tay vào việc vào mùa thu năm ấy, khi các số liệu quan sát của Piazzi được công bố. Chỉ trong vài tuần, ông đã có một quỹ đạo.

Hai ý tưởng

Lời giải của Gauss có hai phần. Phần thứ nhất là một phương pháp mới để xác định quỹ đạo của một thiên thể chỉ từ ba lần quan sát, dựa trên các định luật Kepler: mọi hành tinh chạy trên một hình elip nhận Mặt Trời làm tiêu điểm, theo những quy tắc chính xác về tốc độ. Phần này là một kỳ công tính toán, với rất nhiều phép biến đổi lượng giác, nhưng nó không phải là phần nổi tiếng nhất.

Phần thứ hai là câu hỏi mà mọi nhà khoa học làm việc với số đo đều phải đối mặt: có tới 24 lần quan sát, nhưng lần nào cũng có sai số. Kính thiên văn rung một chút, mắt người đọc vạch chia lệch một chút, không khí làm ảnh sao nhòe đi một chút. Không có quỹ đạo nào đi qua chính xác cả 24 điểm. Vậy nên chọn quỹ đạo nào?

Chỉ dùng ba điểm thì phí mất 21 lần quan sát còn lại, và nếu một trong ba điểm ấy tình cờ lệch nhiều, toàn bộ dự đoán sẽ lệch theo. Lấy trung bình một cách đơn giản thì không biết lấy trung bình cái gì, vì các quan sát ở những thời điểm khác nhau, những vị trí khác nhau.

Gauss trả lời bằng một nguyên tắc: hãy chọn quỹ đạo sao cho tổng bình phương các độ lệch giữa quỹ đạo và các quan sát là nhỏ nhất. Nếu quan sát thứ $i$ lệch khỏi quỹ đạo một khoảng $e_i$ , ta tìm quỹ đạo làm cho

e_1^2 + e_2^2 + \dots + e_n^2

nhỏ nhất có thể. Nguyên tắc này được gọi là phương pháp bình phương nhỏ nhất (method of least squares).

Vì sao lại là bình phương?

Tại sao lại bình phương độ lệch, mà không cộng thẳng độ lệch, hay lấy giá trị tuyệt đối?

Lý do thứ nhất là dấu. Độ lệch có thể âm hoặc dương: quan sát nằm phía trong hay phía ngoài quỹ đạo. Cộng thẳng thì một lệch +5 và một lệch −5 triệt tiêu nhau, cho tổng bằng 0, như thể quỹ đạo khớp hoàn hảo. Bình phương làm mọi độ lệch thành số dương.

Lý do thứ hai là trừng phạt lỗi lớn. Lệch 2 đơn vị bị tính thành 4, lệch 10 đơn vị bị tính thành 100. Phương pháp vì thế rất ghét những quỹ đạo để vài điểm lệch thật xa, và ưu tiên những quỹ đạo lệch một chút với tất cả.

Lý do thứ ba là dễ tính. Đây là lý do thực tế nhất. Với bình phương, bài toán “tìm cực tiểu” biến thành việc giải một hệ phương trình bậc nhất, thứ mà vào năm 1801 có thể làm bằng tay, dù rất lâu.

Gauss còn đi xa hơn. Trong cuốn sách năm 1809, ông chỉ ra rằng nếu sai số đo tuân theo một dạng phân bố hình chuông nhất định, thì bình phương nhỏ nhất cho ra đáp số có khả năng đúng cao nhất[3]. Phân bố hình chuông ấy ngày nay được gọi là phân bố chuẩn, hay phân bố Gauss.

Ceres hiện ra

Các tính toán của Gauss được von Zach công bố cuối năm 1801. Chúng khác hẳn những dự đoán khác. Rồi vào những đêm cuối năm, khi Ceres đã ra khỏi vùng sáng của Mặt Trời, von Zach chĩa kính thiên văn về phía Gauss chỉ, và thấy nó. Ngay sau đó, nhà thiên văn Heinrich Olbers ở Bremen cũng quan sát được. Ceres nằm cách vị trí Gauss dự đoán chưa tới một độ[2].

Chàng trai 24 tuổi lập tức nổi tiếng khắp châu Âu. Năm 1807, Gauss trở thành giám đốc đài thiên văn Göttingen, và giữ chức vụ ấy gần năm mươi năm, tới cuối đời.

Ai nghĩ ra trước?

Gauss không công bố phương pháp của mình ngay. Ông có thói quen chỉ in những gì đã hoàn thiện tới mức không thể bắt bẻ. Phương châm của ông là “ít, nhưng chín”.

Năm 1805, nhà toán học Pháp Adrien-Marie Legendre in một cuốn sách về quỹ đạo sao chổi, trong đó có một phụ lục ngắn trình bày rõ ràng phương pháp bình phương nhỏ nhất, kèm cái tên mà ta dùng đến hôm nay[4]. Bốn năm sau, trong cuốn Theoria motus năm 1809, Gauss trình bày cùng phương pháp và viết rằng ông đã dùng nó từ năm 1795.

Legendre bực bội. Ông viết thư cho Gauss rằng không nên tuyên bố đi trước người khác khi không có ấn phẩm nào chứng minh, và về sau còn công khai phản đối. Gauss chưa bao giờ đưa ra được bằng chứng giấy trắng mực đen. Chỉ có vài đồng nghiệp, trong đó có Olbers, xác nhận rằng Gauss đã cho họ xem phương pháp từ trước năm 1805. Các nhà sử học thống kê ngày nay nhìn chung tin rằng Gauss đã dùng ý tưởng này sớm. Nhưng Legendre là người đầu tiên công bố nó một cách rõ ràng, và Gauss là người đầu tiên đặt nó trên một nền tảng lý thuyết về xác suất[5].

Một mô phỏng nhỏ: vì sao nhiều quan sát thì tốt hơn

Để thấy giá trị của ý tưởng ấy, hãy dựng một mô phỏng đơn giản hóa. Một tiểu hành tinh chạy trên quỹ đạo tròn quanh Mặt Trời, bán kính 2,77 đơn vị thiên văn (AU, khoảng cách Trái Đất – Mặt Trời), giống Ceres. Ta quan sát nó 24 lần trong 41 ngày, mỗi lần đo đều lệch một chút ngẫu nhiên. Rồi ta dự đoán vị trí của nó 300 ngày sau.

Hai cách dự đoán:

Chỉ dùng 3 quan sát: lần đầu, lần giữa và lần cuối.
Dùng cả 24 quan sát và chọn quỹ đạo theo nguyên tắc bình phương nhỏ nhất.

Mỗi cách được thử lại 500 lần, mỗi lần một bộ sai số đo mới:

Cách dự đoán	Sai lệch trung bình	95% trường hợp sai lệch dưới
3 quan sát	0,28 AU	0,66 AU
24 quan sát, bình phương nhỏ nhất	0,13 AU	0,32 AU

Dùng cả 24 quan sát theo cách của Gauss, sai lệch giảm còn khoảng một nửa. Trên bầu trời, đó là khác biệt giữa việc tìm thấy thiên thể ngay trong vùng kính thiên văn đang nhìn, hay phải dò tìm cả một vùng rộng gấp đôi.

Hình tròn nét đứt là quỹ đạo quanh Mặt Trời. Bên phải là cụm 24 chấm xanh, các lần quan sát trong 41 ngày. Phía trên là ngôi sao trắng đánh dấu vị trí thật 300 ngày sau, quanh đó là hai đám chấm dự đoán: đám màu đỏ nhạt từ 3 quan sát trải dài dọc quỹ đạo, đám màu vàng từ bình phương nhỏ nhất gọn hơn hẳn quanh ngôi sao. — Hình 1. Mô phỏng đơn giản hóa (quỹ đạo tròn, đo trực tiếp vị trí), lặp lại 500 lần. Mỗi chấm là một lần dự đoán vị trí sau 300 ngày. Đám màu vàng (bình phương nhỏ nhất với 24 quan sát) tụ sát vị trí thật hơn đám đỏ nhạt (chỉ 3 quan sát). Đây là minh họa ý tưởng, không phải tính toán thật của Gauss. Code ở mục Tự tay tìm lại tiểu hành tinh, gần cuối bài.

Mô phỏng này đơn giản hơn bài toán thật rất nhiều. Gauss phải làm việc với quỹ đạo elip, và với những quan sát chỉ cho biết hướng nhìn từ Trái Đất, cũng đang chuyển động, chứ không cho biết khoảng cách. Nhưng ý tưởng cốt lõi là một: đừng tin vào một vài con số, hãy để mọi quan sát cùng “bỏ phiếu”, và để những sai số ngẫu nhiên tự triệt tiêu lẫn nhau.

Ý tưởng ấy ở đâu hôm nay?

Chấm xanh trên bản đồ

Khi bạn mở bản đồ trên điện thoại, chấm xanh chỉ vị trí của bạn được tính từ tín hiệu của các vệ tinh định vị. Mỗi vệ tinh cho biết khoảng cách từ nó tới điện thoại, nhưng khoảng cách ấy luôn có sai số: tín hiệu bị bầu khí quyển làm chậm, bị nhà cao tầng phản xạ. Về lý thuyết, bốn vệ tinh là đủ để xác định vị trí và thời gian. Nhưng điện thoại thường “nhìn” thấy nhiều vệ tinh hơn thế. Nó không vứt bớt, mà tìm vị trí làm cho tổng bình phương các độ lệch là nhỏ nhất, đúng tinh thần của Gauss với 24 lần quan sát Ceres.

Canh chừng những tiểu hành tinh gần Trái Đất

Bài toán của năm 1801 vẫn đang được giải mỗi ngày, ở quy mô lớn hơn nhiều. Các trung tâm theo dõi tiểu hành tinh tính quỹ đạo của hàng chục nghìn thiên thể bay gần Trái Đất, từ những lần quan sát gửi về từ khắp thế giới, bằng các phiên bản hiện đại của bình phương nhỏ nhất.

Câu chuyện của tiểu hành tinh Apophis cho thấy điều đó rõ nhất. Khi được phát hiện năm 2004, nó chỉ có ít quan sát, và quỹ đạo tính được còn rất bất định. Trong vài ngày cuối năm ấy, khả năng nó va vào Trái Đất năm 2029 được ước tính lên tới vài phần trăm, một con số chưa từng có. Rồi các nhà thiên văn tìm thêm được những bức ảnh chụp nó từ trước, mà trước đó không ai để ý. Có thêm quan sát, vùng bất định co lại, và khả năng va chạm năm 2029 bị loại bỏ hoàn toàn. Năm 2029, Apophis sẽ bay qua rất gần nhưng an toàn. Đúng như mô phỏng ở trên: càng nhiều quan sát, đám mây dự đoán càng gọn.

Đường thẳng khớp nhất

Bình phương nhỏ nhất còn có một tên gọi quen thuộc hơn với sinh viên: hồi quy tuyến tính. Khi ai đó vẽ một đường thẳng “đi qua giữa” một đám điểm dữ liệu, chẳng hạn chiều cao theo cân nặng, hay lượng điện tiêu thụ theo nhiệt độ ngoài trời, gần như chắc chắn họ đang dùng cách của Gauss và Legendre.

Nhiều mô hình học máy (machine learning) hiện đại cũng được huấn luyện theo cùng một tinh thần: chọn tham số sao cho tổng bình phương sai số giữa dự đoán và dữ liệu là nhỏ nhất. Thuật ngữ “hàm mất mát bình phương” (squared loss) mà người làm AI dùng hằng ngày chính là tổng bình phương độ lệch của năm 1801, chỉ là bây giờ có hàng triệu tham số thay vì vài con số mô tả một quỹ đạo.

Tự tay tìm lại tiểu hành tinh

Đoạn Python dưới đây chạy mô phỏng ở trên và vẽ ra Hình 1. Cần cài numpy và matplotlib. Mô hình được đơn giản hóa có chủ ý: quỹ đạo tròn với Mặt Trời ở tâm, và ta đo được trực tiếp vị trí. Bài toán thật của Gauss khó hơn nhiều.

import matplotlib.pyplot as plt
import numpy as np

plt.rcParams.update({"figure.facecolor": "#16181c", "axes.facecolor": "#16181c", "text.color": "#d4d7dc",
                     "axes.edgecolor": "#2a2e35", "xtick.color": "#878d97", "ytick.color": "#878d97", "font.size": 12})
rng = np.random.default_rng(1801)  # Cố định seed để chạy lại ra cùng kết quả.

# Mô phỏng đơn giản hoá: quỹ đạo tròn quanh Mặt Trời, bán kính 2,77 AU, chu kỳ 4,6 năm.
R, PERIOD = 2.77, 4.6 * 365.25
OMEGA = 2 * np.pi / PERIOD                 # tốc độ góc thật, radian mỗi ngày
days = np.linspace(0, 41, 24)              # 24 lần quan sát trong 41 ngày
LATER = 41 + 300                           # cần tìm lại vật thể 300 ngày sau lần quan sát cuối
NOISE = 0.03                               # mỗi lần đo lệch vài phần trăm đơn vị thiên văn


def observe():
    """Vị trí đo được (x, y): vị trí thật cộng sai số đo ngẫu nhiên."""
    angle = OMEGA * days
    xy = np.column_stack([R * np.cos(angle), R * np.sin(angle)])
    return xy + rng.normal(0, NOISE, xy.shape)


def fit(xy, idx):
    """Bình phương nhỏ nhất trên các lần quan sát idx: bán kính r, góc ban đầu, tốc độ góc."""
    r = np.hypot(xy[idx, 0], xy[idx, 1]).mean()
    ang = np.arctan2(xy[idx, 1], xy[idx, 0])
    speed, start = np.polyfit(days[idx], ang, 1)   # đường thẳng khớp nhất với góc theo thời gian
    return r, start, speed


def predict(params, day):
    r, start, speed = params
    return np.array([r * np.cos(start + speed * day), r * np.sin(start + speed * day)])


truth = predict((R, 0.0, OMEGA), LATER)
three, all24 = [], []
for _ in range(500):                               # lặp lại 500 lần, mỗi lần một bộ sai số đo mới
    xy = observe()
    three.append(predict(fit(xy, [0, 12, 23]), LATER))     # chỉ dùng 3 lần quan sát
    all24.append(predict(fit(xy, list(range(24))), LATER)) # dùng cả 24 lần quan sát
three, all24 = np.array(three), np.array(all24)
for name, p in [("3 quan sát", three), ("24 quan sát", all24)]:
    err = np.linalg.norm(p - truth, axis=1)
    print(f"{name:<12} sai lệch trung bình {err.mean():.3f} AU, 95% trường hợp dưới {np.quantile(err, 0.95):.3f} AU")

# Vẽ: quỹ đạo, cung quan sát, và hai "đám mây" dự đoán quanh vị trí thật.
fig, ax = plt.subplots(figsize=(7.5, 7))
t = np.linspace(0, 2 * np.pi, 400)
ax.plot(R * np.cos(t), R * np.sin(t), color="#878d97", lw=1, ls="--", label="quỹ đạo thật")
xy = observe()
ax.plot(xy[:, 0], xy[:, 1], ".", color="#8fa9c9", ms=7, label="24 lần quan sát trong 41 ngày")
ax.plot(three[:, 0], three[:, 1], ".", color="#c98f8f", ms=4, alpha=0.6, label="dự đoán từ 3 quan sát")
ax.plot(all24[:, 0], all24[:, 1], ".", color="#e8b05a", ms=4, alpha=0.8, label="dự đoán bằng bình phương nhỏ nhất (24)")
ax.plot(*truth, "*", color="#d4d7dc", ms=18, label="vị trí thật sau 300 ngày")
ax.plot(0, 0, "o", color="#e8b05a", ms=14)
ax.text(0.15, -0.3, "Mặt Trời", color="#a4a9b2", fontsize=11)
ax.set_aspect("equal")
ax.set_xlim(-3.6, 3.6)
ax.set_ylim(-3.6, 3.6)
ax.legend(loc="center", bbox_to_anchor=(0.5, 0.3), frameon=False, labelcolor="#a4a9b2", fontsize=10)
ax.set_title("Tìm lại một tiểu hành tinh sau 300 ngày (mô phỏng, 500 lần)", color="#d4d7dc", fontsize=12)
fig.savefig("ceres.png", dpi=150, bbox_inches="tight")

3 quan sát   sai lệch trung bình 0.276 AU, 95% trường hợp dưới 0.658 AU
24 quan sát  sai lệch trung bình 0.133 AU, 95% trường hợp dưới 0.319 AU

Hàm np.polyfit chính là bình phương nhỏ nhất: nó tìm đường thẳng (góc theo thời gian) có tổng bình phương độ lệch nhỏ nhất với các quan sát. Bạn có thể thử tăng NOISE hoặc rút ngắn thời gian quan sát để thấy bài toán của Gauss khó tới mức nào khi cung quan sát quá ngắn.

Quay lại đêm cuối năm 1801

Ceres không bị lạc mất. Nó được tìm lại vì một người không chịu vứt đi phần lớn dữ liệu, và cũng không tin mù quáng vào bất kỳ con số đơn lẻ nào. Gauss coi mỗi lần quan sát là một lời chứng có thể sai một chút, rồi tìm câu chuyện phù hợp nhất với tất cả các lời chứng cùng lúc.

Hai trăm năm sau, cách nghĩ ấy có mặt ở khắp nơi: trong chấm xanh trên bản đồ, trong những tính toán canh chừng tiểu hành tinh bay gần Trái Đất, và trong hàng triệu mô hình học máy được huấn luyện mỗi ngày. Năm 2015, tàu vũ trụ Dawn của NASA bay vào quỹ đạo quanh Ceres, chụp được bề mặt của nó từ khoảng cách vài trăm cây số. Cái đốm sáng mờ mà Piazzi suýt để mất giờ là một thế giới có núi, có hố va chạm, có những vệt muối sáng lấp lánh.

Tập trước của loạt bài này kể về Euler và bảy cây cầu Königsberg. Tập sau là về Galois, chàng trai 20 tuổi viết vội những ý tưởng toán học của mình trong đêm trước một cuộc đấu súng.