%% 基本二次函数

x0 = randn(1000,1); % 生成一个1x1000的随机数组，元素值来自标准正态分布
xs = gradient_descent(x0, [1.0], @quadratic_gradient);
disp(all(xs(:,end)==0));

steps = 0.1 * ones(1, 50); % 创建一个包含50个0.1的数组
xs = gradient_descent(x0, steps, @quadratic_gradient); % 调用梯度下降函数

% 计算每个点的二次函数值
ys = arrayfun(@(i) quadratic(xs(:, i)), 1:size(xs, 2));
% 绘制误差曲线
Error_plot(ys, 'log');

%% 约束优化

U = randn(1000, 100);
[U, R] = qr(U); % 计算QR分解，取Q矩阵，如此一来，U就是一个100维空间的正交基
b = randn(1000,1); 

x0 = randn(1000,1);
steps = 0.1 * ones(1, 50);
xs = gradient_descent(x0, steps, @quadratic_gradient, @(x)proj1(x,U,b)); 

x_opt = proj1(0,U,b);%最优解：0的投影

% 绘制迭代点的误差曲线
ys = arrayfun(@(i) quadratic(xs(:, i)), 1:size(xs, 2));
Error_plot(ys, 'linear');
hold on;

% 绘制最优解的误差直线
quadratic_opt = quadratic(x_opt) .* ones(1, size(xs,2));
plot(1:size(xs,2), quadratic_opt, 'r', 'LineWidth', 1); % 绘制红色线
h = legend('误差','$$\frac{1}{2}||x_{opt}||^2$$'); % 添加图例
set(h,'Interpreter','latex');
hold off;

% 绘制迭代点趋近于最优解的过程
ys = arrayfun(@(i) norm(xs(:, i)-x_opt)^2, 1:size(xs, 2));
Error_plot(ys, 'log');


%% 最小二乘法

% 定义基本变量
m = 1000;
n = 100;

A = randn(m,n);
x_opt = randn(n,1);
noise = randn(m,1) * 0.1;
b = A * x_opt + noise;

% 定义目标函数和梯度函数
objective = @(x) least_squares(A, b, x);
gradient = @(x) least_squares_gradient(A, b, x);

% 进行仿真
x0 = randn(n,1);
steps = ones(1,100)*0.1;
xs = gradient_descent(x0, steps, gradient);

% 绘制迭代点的误差曲线
ys = arrayfun(@(i) objective(xs(:, i)), 1:size(xs, 2));
Error_plot(ys, 'log');
hold on;

% 绘制噪声直线
noise_level = norm(noise)^2*ones(1,100);
plot(1:length(steps), noise_level,'r');

% 绘制最优解的误差直线
optimal_level = objective(x_opt)*ones(1,100);
plot(1:length(steps), optimal_level,'g');
legend('error', 'noise level', 'optimal');
hold off;

% 绘制迭代点趋近于最优解的过程
ys = arrayfun(@(i) norm(xs(:, i)-x_opt)^2, 1:size(xs, 2));
Error_plot(ys, 'log');

%% 最小二乘法（m<n）

% 定义基本变量
m = 100;
n = 1000;
A = randn(m,n);
b = randn(m,1);
x_opt = pinv(A)*b;

% 定义目标函数和梯度函数
objective = @(x) least_squares(A, b, x);
gradient = @(x) least_squares_gradient(A, b, x);

% 进行仿真
x0 = randn(n,1);
steps = ones(1,100) * 0.1;
xs = gradient_descent(x0, steps, gradient);

% 绘制迭代点的误差曲线
ys = arrayfun(@(i) objective(xs(:, i)), 1:size(xs, 2));
Error_plot(ys,'log');
hold on;

% 绘制最优解的误差直线
optimal_level=objective(x_opt)*ones(1,100);
plot(1:length(steps),optimal_level,'r');
legend('error', 'optimal');
hold off;

% 绘制迭代点趋近于最优解的过程
ys = arrayfun(@(i) norm(xs(:, i)-x_opt)^2, 1:size(xs, 2));
Error_plot(ys,'linear');
hold on;
optimal_norm=norm(x_opt)^2*ones(1,100);
plot(1:length(steps),optimal_norm,'r');
h=legend('error','$$||x_{opt}||^2$$');
set(h,'Interpreter','latex');

%% l2-正规化

% 定义基本变量
m = 100;
n = 1000;
A = randn(m,n);
b = A*randn(n,1);

% 定义目标函数和梯度函数
objective = @(x) least_squares_l2(A, b, x);
gradient = @(x) least_squares_l2_gradient(A, b, x);

% 理论计算最优解
x_opt = inv(A' * A + 0.1 * eye(1000)) * A' * b;

% 进行仿真
x0 = randn(n,1);
steps = ones(1,500) * 0.1;
xs = gradient_descent(x0, steps, gradient);

% 绘制迭代点的误差曲线
ys = arrayfun(@(i) objective(xs(:, i)), 1:size(xs, 2));
Error_plot(ys,'log');
hold on;

% 绘制最优解的误差直线
optimal_level = objective(x_opt)*ones(1,500);
plot(1:length(steps), optimal_level,'g');
legend('error', 'optimal');
hold off;

% 绘制迭代点趋近于最优解的过程
ys = arrayfun(@(i) norm(xs(:, i)-x_opt)^2, 1:size(xs, 2));
Error_plot(ys,'log');
hold on;
optimal_norm=norm(x_opt)^2*ones(1,500);
plot(1:length(steps),optimal_norm,'r');
h=legend('error','$$||x_{opt}||^2$$');
set(h,'Interpreter','latex');

%% 隐式正则化的魔力

x0 = zeros(n,1);
gradient = @(x) least_squares_gradient(A, b, x);
steps = ones(1,50)*0.1;
xs = gradient_descent(x0,steps,gradient);
ys = arrayfun(@(i) norm(xs(:, i)-x_opt)^2, 1:size(xs, 2));
Error_plot(ys,'linear');
hold on;
optimal_norm=norm(x_opt)^2*ones(1,50);
plot(1:length(steps),optimal_norm,'r');
h=legend('error','$$||x_{opt}||^2$$');
set(h,'Interpreter','latex');

%% LASSO
% 定义基本变量
m = 100;
n = 1000;
A = randn(m,n);
x_opt=zeros(n,1);
x_opt(1:10)=1; % 定义稀疏解
b = A * x_opt;

% 定义目标函数和梯度函数
objective = @(x) lasso(A, b, x);
gradient =@(x) lasso_subgradient(A, b, x);

% 进行仿真
x0 = randn(n, 1);
steps = 0.1 * ones(1, 500); 
xs = gradient_descent(x0, steps, gradient); 

% 绘制迭代点的误差曲线
ys = arrayfun(@(i) objective(xs(:, i)), 1:size(xs, 2));
Error_plot(ys,'log');

% 初始点、最佳点和计算点的比较
idxs = 1:100; % 画前100个点
figure;
plot(idxs, x0(idxs), '--', 'Color', [0.6667 0.6667 0.6667], 'LineWidth', 1.5, 'DisplayName', 'initial');
hold on;
plot(idxs, x_opt(idxs), 'r-', 'LineWidth', 1.5, 'DisplayName', 'optimal');
plot(idxs, xs(end, idxs), 'g-', 'LineWidth', 1.5, 'DisplayName', 'final'); 
xlabel('Coordinate');
ylabel('Value');

%% 支持向量机

z = linspace(-2, 2, 100);
figure('Position', [100, 100, 480, 400]); % 12x10 inches
plot(z, hinge_loss(z), 'LineWidth', 2); 


xlabel('z');
ylabel('Hinge Loss');
title('Hinge Loss Function');

% 定义z值，从-2到2，共100个点
z = linspace(-2, 2, 100);
figure('Position', [100, 100, 480, 400]);
plot(z, hinge_subgradient(z), 'LineWidth', 2); 

xlabel('z');
ylabel('Hinge Subgradient');
title('Hinge Subgradient Function');

% 定义基本变量
m = 1000;
n = 100;
    % 生成m*n的矩阵A，其中前m/2行来自均值为0.1，标准差为1的正态分布，
    % 后m/2行来自均值为-0.1，标准差为1的正态分布
A = [randn(m/2, n)*0.1+randn(m/2, n); -0.1*randn(m/2, n)+randn(m/2, n)];
    % 生成m维向量y，前m/2个元素为1，后m/2个元素为-1
y = [ones(m/2, 1); -1*ones(m/2, 1)];
    % 定义初始点x0
x0 = randn(n, 1);

% 定义目标函数和梯度函数
objective = @(x) svm_objective(A, y, x, 0.05);
gradient = @(x) svm_subgradient(A, y, x, 0.05);

% 执行梯度下降
steps = 0.01 * ones(1, 100);
xs = gradient_descent(x0, steps, gradient); 

% 计算SVM目标函数值的误差曲线
ys = arrayfun(@(i) objective(xs(:, i)), 1:size(xs, 2));

% 绘制误差曲线
Error_plot(ys, 'linear'); 

xavg=0;
for i = 1:size(xs,2)
    xavg = xavg+xs(:,i);
end
xavg = xavg / size(xs,2);
disp("Let's see if averaging out the solutions gives us an improved function value.");
disp(objective(xs(:,end)));
disp(objective(xavg));
    
figure('Position', [100, 100, 480, 400]); % 12x10 inches

ylabel('Accuracy');
xlabel('Step');

% 绘制准确率曲线
hold on; 
plot(1:length(xs), arrayfun(@(i) accuracy(A, y, xs(:,i)), 1:size(xs,2)), 'LineWidth', 2);

% 绘制总体最优解的准确率曲线
plot(1:length(xs), arrayfun(@(i) accuracy(A, y, ones(n, 1)),1:size(xs,2)), 'r', 'LineWidth', 2, 'DisplayName', 'Population optimum');

% 添加图例
legend('Accuracy Curve', 'Population optimum');

%% 稀疏反协方差估计

% 定义基本变量
n = 5;
A = randn(n,n);
S = A*A';

% 定义目标函数
syms Ssyms 5;
syms Xsyms 5;
sparse_inv_cov_syms = trace(Ssyms' * Xsyms) - log(det(Xsyms)) + 0.1 * sum(abs(Xsyms),'all'); % 目标函数
X_vec=Xsyms(:);% 将矩阵 X 向量化
sparse_inv_cov_grad_syms = gradient(sparse_inv_cov_syms,X_vec); % 微分

objective=matlabFunction(sparse_inv_cov_syms,"Vars",{Ssyms,Xsyms}); % 将符号函数转换为matlab函数
grad=matlabFunction(sparse_inv_cov_grad_syms,"Vars",{Ssyms,X_vec});
gradien=@(X) grad(S,X(:));

A0 = randn(n,n);
X0 = A0*A0';
steps=ones(1,500)*0.01;
Xs = gradient_descent_xx(n,X0,steps,gradien,@projection);

ys = arrayfun(@(i) objective(S,Xs(:,(i-1)*n+1:i*n)), 1:500);
Error_plot(ys,'log');


%% 函数定义
function Error_plot(ys, yscale)
% 设置图形大小
    figure('Position', [100, 100, 640, 640]); % 8x8英寸

    % 绘制图形
    plot(1:length(ys), ys,'LineWidth', 2);
    
    % 设置纵坐标的刻度类型
    if strcmpi(yscale, 'log')
        set(gca, 'YScale', 'log');
    else
        set(gca, 'YScale', 'linear');
    end
    
    % 设置横纵坐标标签
    xlabel('Step');
    ylabel('Error');

end

function xs = gradient_descent(init, steps, grad, proj)
    % 投影梯度下降法
    % 输入参数：
    %   init：列向量，迭代起始点
    %   step：标量列表：迭代每一步的步长
    %   grad：函数句柄，梯度函数
    %   proj：函数句柄，投影函数
    % 输出参数：
    %   xs：n行k列矩阵，其中n是init的维数，k是step的长度，每一列表示迭代的每一步
    if nargin < 4
        proj = @(x) x; % proj默认为x
    end
    xs = [init];
    for step = steps
        k=proj(xs(:,end) - step * grad(xs(:,end)));
        xs = [xs, k];
    end
end

function xs = gradient_descent_xx(n,init, steps, grad, proj)
    if nargin < 5
        proj = @(x) x; % If proj is not provided, use the identity function
    end
    xs = [init];
    for step = steps
        startCol = (size(xs, 2) - n) + 1; % 从总列数减去n，然后加1
        endCol = size(xs ,2); % 总列数
        lastSubMatrix = xs(:, startCol:endCol);
        G = grad(lastSubMatrix);
        G = reshape(G,size(lastSubMatrix));
        temp = lastSubMatrix - step * G;
        k=proj(temp);
        xs = [xs, k];
    end
end

function y = quadratic(x)
    y = 0.5 * (x' *x);
end

function grad = quadratic_gradient(x)
    grad = x;
end

function x_proj = proj1(x,U,b)
    x_proj = b + U * (U' * (x-b));
end

function y = least_squares(A, b, x)
    m=size(A,1);
    y=(0.5/m)*norm(A*x-b)^2;
end

function grad = least_squares_gradient(A, b, x)
    m = size(A,1);
    grad = A' * (A*x-b)/m;
end

function obj = least_squares_l2(A, b, x, alpha)
    if nargin < 4
        alpha = 0.1; % 如果没有提供alpha，默认为0.1
    end
    % 最小二乘L2正则化目标函数
    obj = least_squares(A, b, x) + (alpha/2) * x' * x;
end

function grad = least_squares_l2_gradient(A, b, x, alpha)
    if nargin < 4
        alpha = 0.1;
    end
    % 最小二乘L2正则化目标函数的梯度
    grad = least_squares_gradient(A, b, x) + alpha * x;
end

function obj = lasso(A, b, x, alpha)
    if nargin < 4
        alpha = 0.1; 
    end
    % Lasso目标函数
    obj = least_squares(A, b, x) + alpha * norm(x, 1);
end

function g = ell1_subgradient(x)
    % 计算L1范数的次梯度
    g = ones(size(x));
    g(x < 0) = -1;
end

function subgrad = lasso_subgradient(A, b, x, alpha)
    if nargin < 4
        alpha = 0.1;
    end
    % Lasso目标函数的次梯度
    subgrad = least_squares_gradient(A, b, x) + alpha * ell1_subgradient(x);
end

function loss = hinge_loss(z)
    % 返回hinge损失
    loss = max(1.-z, 0);
end

function obj = svm_objective(A, y, x, alpha)
    if nargin < 4
        alpha = 0.1; 
    end
    % SVM目标函数
    m = size(A, 1); 
    obj = mean(hinge_loss(diag(y) * (A * x))) + (alpha/2) * x' * x;
end

function g = hinge_subgradient(z)
    % 计算hinge损失的次梯度
    g = zeros(size(z));
    g(z < 1) = -1;
end

function subgrad = svm_subgradient(A, y, x, alpha)
    if nargin < 4
        alpha = 0.1; 
    end
    % 计算SVM目标函数的次梯度
    g1 = hinge_subgradient(diag(y) * (A * x));
    g2 = diag(y) * A;
    subgrad = (g1' * g2)' + alpha * x;
end

function acc = accuracy(A, y, x)
    % 计算准确率
    acc = mean(diag(y) * (A * x) > 0);
end

function obj = sparse_inv_cov(S, X, alpha)
    if nargin < 3
        alpha = 0.1; 
    end
    obj = trace(S' * X) - log(det(X)) + alpha * sum(abs(X));
end

function X_proj = projection(X)
    % Projection onto positive semidefinite cone.
    [U,es] = eig(X);
    es(es < 0) = 0;
    X_proj = U * es * U';
end
