![logo](GCTL_logo.jpg)

# GCTL Optimization Library

> **注意**: 本文档由 Cursor AI 自动生成。文档中的函数名称、参数和用法说明可能存在错误，请以头文件（`.h`）中的实际函数声明为准。如发现任何不一致，请以源代码为准。

GCTL Optimization Library 是一个高效的数值优化算法库，提供了多种优化算法实现和求解器。该库是 Generic Scientific Template Library (GSTL) 的一部分。

## 主要特性

- 多种优化算法实现
  - LCG (Linear Conjugate Gradient): 线性共轭梯度法
  - CLCG (Complex Linear Conjugate Gradient): 复数域线性共轭梯度法
  - L-BFGS (Limited-memory BFGS): 限制内存BFGS算法
  - SGD (Stochastic Gradient Descent): 随机梯度下降
  - LGD (Lévy Gradient Descent): 基于Lévy飞行的随机搜索方法，具有全局收敛性
  - SVD (Singular Value Decomposition): 奇异值分解

- 矩阵分解工具
  - Cholesky分解
  - LU分解
  - SVD分解

- 损失函数支持
  - 提供多种常用损失函数
  - 支持自定义损失函数

## 安装

项目使用CMake构建系统。安装步骤如下：

```bash
mkdir build
cd build
cmake ..
make
```

## 使用方法

### 1. 配置求解器

所有求解器都支持通过TOML配置文件设置参数：

```toml
[solver]
max_iterations=1000    # 最大迭代次数
epsilon=1e-6          # 收敛阈值
abs_diff=1            # 是否使用绝对差异
restart_epsilon=1e-4  # 重启阈值
step=0.1             # 步长
sigma=0.9            # sigma参数
beta=0.1             # beta参数
maxi_m=10            # 最大M值（用于L-BFGS）
```

### 2. 定义优化问题

GCTL优化库通过类的虚函数接口来定义具体的优化问题。不同类型的优化问题需要实现不同的虚函数接口：

#### 2.1 线性问题接口

对于求解线性方程组 Ax = b 的问题：

```cpp
class LinearProblem : public gctl::lcg_solver {
public:
    // 必须实现：定义矩阵向量乘积 Ax
    virtual void LCG_Ax(const gctl::array<double> &x, 
                       gctl::array<double> &ax) override {
        // 示例：实现 ax = A * x
        for(size_t i = 0; i < n; ++i) {
            ax[i] = 0.0;
            for(size_t j = 0; j < n; ++j) {
                ax[i] += A[i][j] * x[j];
            }
        }
    }
    
    // 可选实现：定义预处理矩阵M的作用
    virtual void LCG_Mx(const gctl::array<double> &x, 
                       gctl::array<double> &mx) override {
        // 示例：对角预处理
        for(size_t i = 0; i < n; ++i) {
            mx[i] = x[i] / A[i][i];
        }
    }
};
```

#### 2.2 非线性优化问题接口

对于求解非线性优化问题 min f(x) 的情况：

```cpp
class NonlinearProblem : public gctl::lbfgs_solver {
public:
    // 必须实现：计算目标函数值和梯度
    virtual double LBFGS_Evaluate(const gctl::array<double> &x,
                                gctl::array<double> &g) override {
        // 示例：Rosenbrock函数及其梯度
        double fx = 0.0;
        for(size_t i = 0; i < n-1; ++i) {
            fx += 100.0 * pow(x[i+1] - x[i]*x[i], 2) + 
                  pow(1.0 - x[i], 2);
        }
        
        // 计算梯度
        for(size_t i = 0; i < n; ++i) {
            if(i == 0) {
                g[i] = -400.0*x[i]*(x[i+1] - x[i]*x[i]) - 
                       2.0*(1.0 - x[i]);
            }
            else if(i == n-1) {
                g[i] = 200.0*(x[i] - x[i-1]*x[i-1]);
            }
            else {
                g[i] = 200.0*(x[i] - x[i-1]*x[i-1]) - 
                       400.0*x[i]*(x[i+1] - x[i]*x[i]) - 
                       2.0*(1.0 - x[i]);
            }
        }
        return fx;
    }
    
    // 可选实现：定义预处理
    virtual void LBFGS_Precondition(const gctl::array<double> &x,
                                  const gctl::array<double> &g,
                                  const gctl::array<double> &d,
                                  gctl::array<double> &d_pre) override {
        // 实现预处理
        d_pre = d;  // 默认不做预处理
    }
};
```

#### 2.3 全局优化问题接口

对于需要进行全局搜索的优化问题：

```cpp
class GlobalProblem : public gctl::lgd_solver {
public:
    // 必须实现：计算目标函数值和梯度
    virtual double LGD_Evaluate(const gctl::array<double> &x,
                              gctl::array<double> &g) override {
        // 示例：多峰函数及其梯度
        double fx = 0.0;
        for(size_t i = 0; i < n; ++i) {
            fx += sin(x[i]) * exp(-0.1 * x[i]*x[i]);
        }
        
        // 计算梯度
        for(size_t i = 0; i < n; ++i) {
            g[i] = cos(x[i]) * exp(-0.1 * x[i]*x[i]) - 
                   0.2 * x[i] * sin(x[i]) * exp(-0.1 * x[i]*x[i]);
        }
        return -fx;  // 最小化负值等价于最大化原函数
    }
    
    // 可选实现：监控优化进度
    virtual int LGD_Progress(const int curr_t, const double curr_fx,
                           const double mean_fx, const double best_fx,
                           const lgd_para &param) override {
        // 返回非零值将终止优化
        return 0;
    }
};
```

#### 2.4 随机优化问题接口

对于需要随机采样或批量处理的优化问题：

```cpp
class StochasticProblem : public gctl::sgd_solver {
public:
    // 必须实现：计算单个批次的损失
    virtual double evaluate_batch(const gctl::array<double> &x,
                                const gctl::array<size_t> &batch_indices) override {
        // 示例：均方误差损失
        double loss = 0.0;
        for(auto idx : batch_indices) {
            double pred = predict(x, data[idx]);
            loss += 0.5 * pow(pred - target[idx], 2);
        }
        return loss / batch_indices.size();
    }
    
    // 必须实现：计算单个批次的梯度
    virtual void gradient_batch(const gctl::array<double> &x,
                              const gctl::array<size_t> &batch_indices,
                              gctl::array<double> &grad) override {
        // 示例：均方误差损失的梯度
        grad.fill(0.0);
        for(auto idx : batch_indices) {
            double pred = predict(x, data[idx]);
            update_gradient(x, data[idx], pred - target[idx], grad);
        }
        grad *= (1.0 / batch_indices.size());
    }
};
```

### 3. 基本用法

#### 2.1 线性共轭梯度法 (LCG)

```cpp
#include <gctl/optimization.h>

// 继承lcg_solver类并实现必要的虚函数
class MyProblem : public gctl::lcg_solver {
public:
    // 实现矩阵向量乘法 Ax
    virtual void LCG_Ax(const gctl::array<double> &x, gctl::array<double> &ax) override {
        // 实现 ax = A * x
    }
    
    // 实现预处理矩阵乘法 Mx（可选）
    virtual void LCG_Mx(const gctl::array<double> &x, gctl::array<double> &mx) override {
        // 实现预处理，如果不需要预处理，可以简单地将mx = x
        mx = x;
    }
};

// 使用求解器
void solve() {
    MyProblem solver;
    solver.set_max_iterations(1000);
    solver.set_epsilon(1e-6);
    solver.solve(b, x);  // b为右端向量，x为解向量
}
```

#### 2.2 复数域线性共轭梯度法 (CLCG)

```cpp
#include <gctl/optimization.h>

class ComplexProblem : public gctl::clcg_solver {
public:
    virtual void CLCG_Ax(const gctl::array<std::complex<double>> &x,
                        gctl::array<std::complex<double>> &ax) override {
        // 实现复数域的矩阵向量乘法
    }
};
```

#### 2.3 L-BFGS优化器

```cpp
#include <gctl/optimization.h>

class OptProblem : public gctl::lbfgs_solver {
public:
    // 计算目标函数值
    virtual double evaluate(const gctl::array<double> &x) override {
        // 返回在点x处的函数值
    }
    
    // 计算梯度
    virtual void gradient(const gctl::array<double> &x, 
                         gctl::array<double> &grad) override {
        // 计算在点x处的梯度
    }
};
```

#### 2.4 Lévy梯度下降 (LGD)

```cpp
#include <gctl/optimization.h>

class LGDProblem : public gctl::lgd_solver {
public:
    // 计算目标函数值
    virtual double evaluate(const gctl::array<double> &x) override {
        // 返回在点x处的函数值
    }
    
    // 计算梯度
    virtual void gradient(const gctl::array<double> &x, 
                         gctl::array<double> &grad) override {
        // 计算在点x处的梯度
    }
    
    // 设置Lévy飞行参数
    void setup() {
        set_levy_alpha(1.5);  // 设置Lévy分布的α参数
        set_step_size(0.1);   // 设置步长
    }
};
```

#### 2.5 随机梯度下降 (SGD)

```cpp
#include <gctl/optimization.h>

class SGDProblem : public gctl::sgd_solver {
public:
    // 计算随机批次的损失
    virtual double evaluate_batch(const gctl::array<double> &x,
                                const gctl::array<size_t> &batch_indices) override {
        // 返回当前批次的损失值
    }
    
    // 计算随机批次的梯度
    virtual void gradient_batch(const gctl::array<double> &x,
                              const gctl::array<size_t> &batch_indices,
                              gctl::array<double> &grad) override {
        // 计算当前批次的梯度
    }
    
    // 配置SGD参数
    void setup() {
        set_batch_size(32);           // 设置批次大小
        set_learning_rate(0.01);      // 设置学习率
        set_momentum(0.9);            // 设置动量（可选）
    }
};
```

### 3. 高级功能

#### 3.1 自定义损失函数

```cpp
#include <gctl/optimization.h>

class CustomLoss : public gctl::loss_function {
public:
    virtual double evaluate(const gctl::array<double> &pred,
                          const gctl::array<double> &target) override {
        // 实现自定义损失函数计算
    }
    
    virtual void gradient(const gctl::array<double> &pred,
                         const gctl::array<double> &target,
                         gctl::array<double> &grad) override {
        // 实现损失函数关于预测值的梯度计算
    }
};
```

#### 3.2 矩阵分解工具使用

```cpp
#include <gctl/optimization.h>

void matrix_operations() {
    // SVD分解
    gctl::matrix<double> A, U, S, V;
    gctl::svd_decomp(A, U, S, V);
    
    // LU分解
    gctl::matrix<double> L, U;
    gctl::lu_decomp(A, L, U);
    
    // Cholesky分解
    gctl::matrix<double> L;
    gctl::cholesky_decomp(A, L);
}
```

更多详细示例请参考 `example` 目录下的示例代码。

## 示例说明

- `ex1.cpp`: 基本的共轭梯度求解示例
- `ex2.cpp`: L-BFGS优化示例
- `ex3.cpp`: SVD分解应用
- `ex4.cpp`: 随机梯度下降示例
- `ex5.cpp`: 约束优化问题
- `ex6.cpp`: 莱维梯度下降示例
- `ex7.cpp-ex10.cpp`: 更多高级应用示例

## 许可证

该项目采用双重许可：
1. GNU Lesser General Public License v2.0
2. 商业许可（需单独申请）

## 联系方式

作者：Yi Zhang (yizhang-geo@zju.edu.cn)

## 参与贡献

欢迎提交Issue和Pull Request来帮助改进项目。