Mrówka na sferze: klasyczne sterowanie R³ kontra optymalizacja riemannowska na S²

Problem #

Dana jest sfera S². Na sferze siedzi mrówka, którą możemy sterować. Chcemy, żeby wykonała następującą trasę:

• Start: punkt p₀ = (1, 0, 0) na równiku
• Etap 1: dojść do bieguna p₁ = (0, 0, 1) najkrótszą drogą
• W biegunie: obrócić się o 90° w miejscowej przestrzeni stycznej T_{p₁}S²
• Etap 2: wrócić na równik do punktu p₂ = (0, 1, 0)

Funkcja kosztu: mrówka ma poruszać się bez szarpania — minimalizujemy całkę z kowariantnego przyspieszenia wzdłuż trajektorii.

Na pozór prosty problem. A jednak podejście klasyczne — optymalizacja w R³ — napotka tu fundamentalne trudności, których żadna kara lagrange’owska nie usunie.

Holonomia: dwie różne niespodzianki #

Warto tu oddzielić dwie rzeczy, które łatwo pomylić.

Co jest zadane wprost #

W naszym przykładzie mrówka aktywnie obraca się o 90° na biegunie — to jest warunek brzegowy który zadaliśmy wprost. Wektor prędkości przy przybyciu do bieguna wynosi (-1, 0, 0); obracamy go do (0, 1, 0) i stamtąd mrówka wyrusza z powrotem. Punkt końcowy (0, 1, 0) leży 90° od punktu startowego (1, 0, 0) na równiku — ale to konsekwencja naszego wyboru obrotu, nie efekt krzywizny.

Co wynika z krzywizny sfery #

Holonomia to inny, choć pokrewny efekt. Wyobraź sobie że mrówka startuje z (1, 0, 0) niosąc wektor wskazujący na biegun — i porusza się po tym samym trójkącie geodezyjnym, nie obracając wektora aktywnie w żadnym punkcie. Wektor jest transportowany równolegle: w każdym punkcie trasy pozostaje “jak najbardziej stały” w sensie lokalnym, bez aktywnego skręcania.

Gdy mrówka wróci na równik do (1, 0, 0), jej wektor będzie obrócony względem pierwotnego. O ile? Dokładnie o kąt równy całce krzywizny Gaussa po obszarze otoczonym przez trasę:

Ω = ∫∫ K dA = K · A = 1 · (π/2) = π/2

gdzie K = 1 to krzywizna Gaussa sfery jednostkowej, a A = π/2 to pole trójkąta geodezyjnego. Wektor obrócił się o 90° — choć nikt go nie kręcił. To jest właśnie holonomia sfery — efekt globalny wynikający z lokalnej krzywizny, niewidoczny żadnemu mieszkańcowi sfery poruszającemu się „po prostej".

W naszym przykładzie z mrówką te dwa efekty są splecione: aktywny obrót w biegunie jest zadany, a krzywizna sfery wpływa na geometrię geodezyjnych i strukturę przestrzeni stycznej wzdłuż trasy. Rozdzielenie ich jest ważne, bo holonomia działa nawet bez żadnego aktywnego obrotu — wystarczy zamknięta krzywa na zakrzywionej powierzchni.

Podejście 1: optymalizacja w R³ #

Sformułowanie #

Dyskretyzujemy trajektorię jako N punktów {x₀, x₁, …, x_{N-1}} ⊂ R³. Minimalizujemy dyskretne przyspieszenie:

C_R³ = Σᵢ ‖xᵢ₊₁ - 2xᵢ + xᵢ₋₁‖²

Więz ‖xᵢ‖ = 1 dodajemy jako karę:

C = C_R³ + λ · Σᵢ (‖xᵢ‖² - 1)²

Optymalizacja przez scipy.optimize.minimize (L-BFGS-B).

Problemy, które napotkamy #

Problem 1 — inicjalizacja wychodzi ze sfery. Liniowa interpolacja między dwoma punktami na sferze daje prostą w R³, nie krzywą na sferze. Punkt środkowy między (1,0,0) a (0,0,1) to (0.5, 0, 0.5) z normą ≈ 0.71. Inicjalizacja jest już „pod ziemią".

Problem 2 — optymalizator bez więzu osiąga koszt zero. Minimalizacja C_R³ bez kary daje wynik idealny — ale trajektoria zapada się do środka sfery. Najgładszą krzywą w R³ między dwoma punktami jest odcinek, i właśnie tak optymalizator rozwiązuje problem. Koszt = 0.0, mrówka pod ziemią.

Problem 3 — paradoks geodezyjnych. To jest sedno sprawy. Geodezyjna na sferze to wielkie koło — jest zakrzywiona w R³, bo podąża za krzywizną powierzchni. Jej dyskretne przyspieszenie euklidesowe jest niezerowe:

Koszt C_R³ geodezyjnej:      0.075540
Koszt C_R³ opt. bez więzu:   0.000000

Metryka R³ karze geodezyjne i nagradza zejście ze sfery. Cel optymalizacji i geometria problemu są fundamentalnie sprzeczne.

Problem 4 — kara jako kompromis. Przy λ = 10 normy punktów wychodzą 0.9996–1.0000, błąd więzu ≈ 4·10⁻⁴. Niby nieźle, ale to wynik kompromisu między dwiema walczącymi ze sobą funkcjami: „bądź gładki w R³" vs „zostań na sferze". Żadna z nich nie rozumie geometrii S².

Wyniki numeryczne (R³) #

Podejście 2: optymalizacja riemannowska na S² #

Przestrzeń konfiguracji #

Zamiast optymalizować punkty w R³ z karą, wprost pracujemy na produkcie kartezjańskim rozmaitości:

M = S² × S² × ... × S²   (M = N-2 razy, tylko punkty wewnętrzne)

Pymanopt obsługuje tę strukturę przez Product([Sphere(3)] * M).

Funkcja kosztu: przyspieszenie kowariantne #

Zamiast różnic skończonych w R³, używamy odwzorowania logarytmicznego na sferze:

Log_x(y) = arccos(⟨x,y⟩) · (y - ⟨x,y⟩·x) / ‖y - ⟨x,y⟩·x‖

Jest to wektor w T_x S² wskazujący na y — „prędkość", z jaką geodezyjna startująca z x dociera do y. Dyskretne przyspieszenie kowariantne w punkcie xᵢ:

∇_γ̇ γ̇ ≈ Log_{xᵢ}(xᵢ₊₁) + Log_{xᵢ}(xᵢ₋₁)

To jest dyskretna wersja drugiej pochodnej kowariantnej — miara „szarpania" mierzona na sferze, nie w przestrzeni otaczającej. Geodezyjna ma ∇_γ̇ γ̇ = 0 wszędzie, bo jest „krzywą bez przyspieszenia" w sensie riemannowskim.

Funkcja kosztu:

@pymanopt.function.autograd(manifold)
def cost(*args):
    x = build_full(np.array(args))
    total = 0.0
    for i in range(1, N - 1):
        v_fwd = log_map(x[i], x[i+1])
        v_bwd = log_map(x[i], x[i-1])
        acc_cov = v_fwd + v_bwd
        total += np.dot(acc_cov, acc_cov)
    return total

Backend autograd oblicza gradient euklidesowy automatycznie. Pymanopt rzutuje go na T_x S² i wykonuje retrakcję — cała mechanika riemannowska jest przezroczysta dla użytkownika.

Co dzieje się inaczej niż w R³ #

Inicjalizacja: geodezyjne na S². Normy wszystkich punktów: 1.00000000 / 1.00000000. Nie ma co rzutować — punkty leżą na sferze z definicji.

Optymalizacja: gradient zawsze leży w T_x S², retrakcja nigdy nie opuszcza sfery. Więz nie istnieje jako ograniczenie — jest wbudowany w geometrię przestrzeni konfiguracji.

Błąd więzu po optymalizacji: 1.11·10⁻¹⁶ — epsilon maszynowy. Nie jest to kompromis; to dokładność arytmetyki zmiennoprzecinkowej.

Wyniki numeryczne (S²) #

Koszt zmalał o ponad dwa rzędy wielkości. Pozostały niezerowy koszt wynika z kąta między dwiema geodezyjnymi w biegunie — to jest geometryczna osobliwość problemu, nie błąd algorytmu. Optymalizator wygładził trajektorię tak bardzo, jak to możliwe przy zadanych warunkach brzegowych.

Zestawienie obu podejść #

Dlaczego to ma znaczenie poza matematyką #

Ten przykład jest minimalnym przybliżeniem klasy problemów inżynierskich. Wszędzie tam, gdzie przestrzeń stanów jest naturalnie rozmaitością — orientacje w SO(3), układy dq w maszynach elektrycznych, stany w mechanice kwantowej — podejście R³ z karami wprowadza te same patologie:

• inicjalizacja opuszcza rozmaitość
• gradient nie jest tangentny
• metryka przestrzeni otaczającej nie odpowiada metryce problemu
• wynik zależy od strojenia parametrów kary

Podejście riemannowskie usuwa wszystkie cztery problemy jednocześnie — nie przez dodanie ograniczeń, lecz przez zmianę przestrzeni, w której pracujemy.

Podziękowania #

Przykład korzysta z biblioteki Pymanopt — narzędzia do optymalizacji na rozmaitościach riemannowskich z automatycznym różniczkowaniem. Jeśli korzystasz z pymanopt w pracy naukowej, cytuj:

James Townsend, Niklas Koep, Sebastian Weichwald. Pymanopt: A Python Toolbox for Optimization on Manifolds using Automatic Differentiation. Journal of Machine Learning Research, 17(137):1–5, 2016. http://jmlr.org/papers/v17/16-177.html

@article{JMLR:v17:16-177,
    author  = {James Townsend and Niklas Koep and Sebastian Weichwald},
    journal = {Journal of Machine Learning Research},
    number  = {137},
    pages   = {1--5},
    title   = {Pymanopt: A Python Toolbox for Optimization on Manifolds
               using Automatic Differentiation},
    url     = {http://jmlr.org/papers/v17/16-177.html},
    volume  = {17},
    year    = {2016}
}