let menu = ["Home", "Algorithms", "CodeHub", "VNOI Statistics"];

Tìm cực tiểu hàm số bằng Gradient descent

Tìm cực tiểu hàm số bằng Gradient descent

Gradient descent là phương pháp dùng để tìm cực tiểu địa phương của hàm số. Nó có thể được áp dụng trong trường hợp đại lượng cần tìm có thể được biểu diễn thành hàm số.

Phần sau đây sẽ trình bày các kiến thức cần thiết để sử dụng Gradient descent.

Đạo hàm riêng

Định nghĩa

Cho \(f\) là hàm số hai biến \(x\) và \(y\). Nếu ta xem \(y\) như hằng số và lấy đạo hàm theo \(x\), ta được đạo hàm riêng của \(f\) theo \(x\), ký hiệu bởi \(f_x\) . Tương tự, ta được đạo hàm riêng của \(f\) theo \(y\) , ký hiệu bởi \(f_y\) . Nói cách khác, \(f_x\) và \(f_y\) là các hàm số được định bởi:

\[f_x(x, y) = \lim_{h \to 0} \frac{f(x+h, y) - f(x, y)}{h}\] \[f_y(x, y) = \lim_{h \to 0} \frac{f(x, y+h) - f(x, y)}{h}\]

Miễn là các giới hạn trên tồn tại.

Các công thức cho đạo hàm một biến vẫn có thể được áp dụng.

Các kí hiệu khác

Nếu viết \(z = f(x, y)\), người ta cũng có nhiều ký hiệu khác cho đạo hàm riêng như sau:

\[f_x = \frac{\delta f}{\delta x} = \frac{\delta z}{\delta x} = f_1 = D_1 f = D_x f\] \[f_y = \frac{\delta f}{\delta y} = \frac{\delta z}{\delta y} = f_2 = D_2 f = D_y f\]

Đạo hàm cho hàm số nhiều biến được định nghĩa và kí hiệu tương tự.

Ví dụ

Với \(f(x, y) = x^3 + x^2y^3 - 2y^2\), xem \(y\) là hằng số, lấy đạo hàm theo \(x\) ta được:

\[f_x(x, y) = 3x^2 + 2xy^3\]

Tương tự ta có:

\[f_y(x, y) = 3x^2y^2 - 4y\]

Gradient

Định nghĩa

Gradient của hàm số \(f(x, y)\) là một trường vector được định nghĩa như sau:

\[\nabla f(x, y) = \langle f_x(x, y), f_y(x, y) \rangle\]

Nói cách khác, gradient của hàm số \(f\) là một hàm vector \(\nabla f\) với mỗi thành phần lần lượt là các đạo hàm theo biến \(x\) và theo biến \(y\) của \(f\).

Tính chất

Xét hàm 2 biến \(f\) tại điểm \((x, y)\), nếu ta đi theo hướng của \(\nabla f\) thì giá trị của hàm số sẽ thay đổi nhanh nhất.

Một cách trực quan hơn, nếu xem đồ thị của hàm số 2 biến \(f\) là một vùng đồi núi thì tại điểm có toạ độ \((x, y)\), hướng \(\nabla f\) sẽ là hướng có độ dốc cao nhất.

Gradient descent

Thuật toán Gradient descent hoạt động dựa trên nhận xét: Nếu hàm số \(f(x)\) xác định và khả vi tại điểm \(a\), khi đó giá trị của \(f\) sẽ giảm nhanh nhất khi đi theo hướng ngược với gradient của \(f\) tại a, \(- \nabla f(a)\). Với:

\[\large a_{n+1} = a_n - \gamma \nabla f(a_n)\]

Và \(\gamma\) đủ nhỏ thì ta sẽ có \(f(a_n) \ge f(a_{n+1})\).

Nói cách khác, nếu ta chọn điểm xuất phát \(x_0\), và sau đó đi theo công thức trên thì ta sẽ đi dần đến điểm cực tiểu (local minimum) của hàm \(f\) (cực tiểu ở đây không chắc là giá trị nhỏ nhất của hàm số).

Đoạn code python sau tìm cực tiểu của hàm \(f(x) = x^4 - 3x^3 + 2\), có đạo hàm là \(f'(x) = 4x^3 - 9x^2\) (đối với hàm số 1 biến, gradient của hàm số chính là đạo hàm):

# From calculation, it is expected that the local minimum occurs at x=9/4

cur_x = 6 # The algorithm starts at x=6
gamma = 0.01 # step size multiplier
precision = 0.00001
previous_step_size = cur_x

def df(x):
    y = 4 * x**3 - 9 * x**2
    return y

while previous_step_size > precision:
    prev_x = cur_x
    cur_x += -gamma * df(prev_x)
    previous_step_size = abs(cur_x - prev_x)

print("The local minimum occurs at %f" % cur_x)

Khi chạy thử đoạn code trên, nếu ta chọn gamma quá lớn, thuật toán sẽ không chính xác. Mặt khác, nếu chọn gamma quá nhỏ thì thuật toán sẽ rất chậm. Để chọn \(\gamma\) cho phù hợp, ta có công thức Barzilai-Borwein:

\[\gamma _n = \frac{(x_n - x_{n-1})^T[\nabla f(x_n) - \nabla f(x_{n-1})]}{\|\nabla f(x_n) - \nabla f(x_{n-1})\|^2}\]

Ứng dụng

Bài toán, phân tích

Cho \(n\) điểm trên mặt phẳng toạ độ \(Oxy\), tìm điểm \(X\) sao cho tổng khoảng cách từ \(X\) đến \(n\) điểm đã cho là nhỏ nhất.

Bài toán trên có tên là Geometric median. Để giải, ta có thể quy về tìm cực tiểu của hàm số 2 biến:

\(f(X) = \displaystyle\sum_{i=1}^{n} \|X - A_i\|\) , với \(A_i\) là những điểm đã cho.

Ta sẽ dùng phương pháp gradient descent để giải bài toán trên. Đầu tiên, cần tìm các đạo hàm riêng của hàm số, ta có:

Biễu diễn lại hàm số trên:

\[f(X) = f(x, y) = \displaystyle \sum \sqrt{(x-x_i)^2 + (y-y_i)^2}\]

Đạo hàm riêng:

\[f_x = \sum \frac{x - x_i}{\sqrt{(x-x_i)^2 + (y-y_i)^2}}\] \[f_y = \sum \frac{y - y_i}{\sqrt{(x-x_i)^2 + (y-y_i)^2}}\]

Gradient:

\[\nabla f(X) = \langle f_x, f_y \rangle = \sum \frac{X - A_i}{\| X - A_i \|}\]

Cài đặt

Khi cài đặt ta sẽ chọn điểm đầu là trọng tâm của các điểm đang xét để thuật toán có thể chạy tốt.

#include <cmath>
#include <vector>
#include <iostream>
#include <sstream>
using namespace std;

struct vec {
    double x, y;
    vec& operator+=(vec v) {
        this->x += v.x, this->y += v.y;
        return *this;
    }
    vec& operator-=(vec v) {
        this->x -= v.x, this->y -= v.y;
        return *this;
    }
    vec& operator*=(double k) {
        this->x *= k, this->y *= k;
        return *this;
    }
    vec& operator/=(double k) {
        this->x /= k, this->y /= k;
        return *this;
    }
};
vec operator+(vec u, vec v) { return u += v; }
vec operator-(vec u, vec v) { return u -= v; }
double operator*(vec u, vec v) { return u.x*v.x + u.y*v.y; }
vec operator*(vec v, double u) { return v *= u; }
vec operator/(vec v, double u) { return v /= u; }
double magnitude(vec u) { return sqrt(u.x*u.x + u.y*u.y); }
double sqr(double d) { return d*d; }

vec gradient(const vector<vec> &a, vec x) {
    vec tmp, res = {0, 0};
    for (auto &u: a) {
        tmp = x - u;
        res += tmp / magnitude(tmp);
    }
    return res;
}

double f(const vector<vec> &a, vec x) {
    double res = 0;
    for (auto& u: a) res += magnitude(x - u);
    return res;
}

double compute(const vector<vec> &a) {
    const double epsilon = 1e-7;
    double gamma = 0.01;
    vec x0, x1 = {0, 0};
    for (auto& u: a) x1 += u;
    x1 /= (double) a.size();
    double f0, f1 = f(a, x1);
    vec g0, g1 = gradient(a, x1);
    do {
        x0 = x1, x1 -= g1 * gamma;
        f0 = f1, f1 = f(a, x1);
        g0 = g1, g1 = gradient(a, x1);
        gamma = (x1 - x0) * (g1 - g0) / sqr(magnitude(g1 - g0));
    } while (abs(f1 - f0) > epsilon);
    return f1;
}

int main() {
    // Input: dòng đầu tiên là số điểm N, N dòng tiếp theo
    // mỗi dòng gồm 2 số x y là toạ độ của điểm.
    ios::sync_with_stdio(false), cin.tie(0);
    int n; cin >> n;
    vector<vec> a(n);
    for (auto& u: a) cin >> u.x >> u.y;
    cout << compute(a) << "\n";
    return 0;
}
Comments