如何按给定算法手动编写Matlab 1范数回归自定义函数

现有代码的明显错误 #

• 第6行语法错误：n = size(1,A); 应修正为 n = size(A,2);，size函数第一个参数为目标矩阵，第二个参数为维度标识，1对应行数、2对应列数，原写法会直接触发运行报错。
• 第20行变量名错误：y_s = y(s) 应修正为 y_s = y_vect(s);，你代码中定义的对偶向量变量名为y_vect，原写法会报变量未定义错误，同时末尾缺少分号会在命令行打印无关输出。
• 无迭代循环结构：你仅写了单次迭代的逻辑，且else分支中计算出的x_new没有赋值给x_temp，代码最终只会返回初始计算的x_temp，完全无法实现迭代优化效果。
• 数值稳定性问题：多处调用inv()求解矩阵逆，计算效率低且容易出现数值误差，建议替换为MATLAB矩阵左除运算符\。
• 初始基集校验缺失：直接默认矩阵前n列为可逆基，没有做可逆性校验，遇到前n列线性相关的输入会直接报奇异矩阵错误。

算法迭代逻辑说明 #

你实现的是求解切比雪夫极小化问题 $\min_{x} |Ax - b|_\infty$ 的交换算法，完整迭代逻辑如下：

1. 初始选取n个样本点构成基集set_B，校验基集对应子矩阵可逆后，求解基集下的方程组得到初始解x_temp
2. 计算所有样本点的残差 h = A*x_temp - b
3. 构造对偶向量y_vect：基外样本点的对偶值取对应残差的符号，基内样本点的对偶值通过子矩阵运算求得
4. 最优性校验：如果基内所有样本点对应的对偶分量绝对值都≤1，当前解就是全局最优解，直接返回
5. 不满足最优条件时，选取第一个对偶分量绝对值超过1的基内索引s，计算步长向量t_vect
6. 基外样本点中筛选出正步长，取最小步长对应的样本索引r，将s移出基集、r加入基集完成基交换
7. 用新的基集重新计算解，回到步骤2循环迭代直到满足最优条件

修正后可运行代码 #

function x_opt = chebyshev_solver(A,b)
% 求解切比雪夫极小化问题 min_x ||Ax - b||_inf
% 输入：A m*n系数矩阵，b m维列向量
% 输出：最优解x_opt n维列向量
    b_vect = b(:); % 统一转为列向量
    [m,n] = size(A);
    if m < n
        error('样本数不能小于变量维数');
    end
    % 初始基集选取：优先选前n列，若不可逆则随机选取直到可逆
    set_B = 1:n;
    while rank(A(set_B,:)) < n
        set_B = randperm(m,n);
    end
    set_B_Comp = setdiff(1:m, set_B);
    max_iter = 1000; % 最大迭代次数防死循环
    iter = 0;
    x_opt = [];
    while iter < max_iter
        iter = iter + 1;
        % 求解当前基下的解
        x_temp = A(set_B, :) \ b_vect(set_B);
        h = A*x_temp - b_vect;
        % 计算对偶向量
        y_vect = zeros(m, 1);
        y_vect(set_B_Comp) = sign(h(set_B_Comp));
        y_vect(set_B) = - (A(set_B, :)') \ (A(set_B_Comp, :)' * y_vect(set_B_Comp));
        % 最优性校验，加容差处理数值误差
        abs_y_B = abs(y_vect(set_B));
        if all(abs_y_B <= 1 + 1e-6) 
            x_opt = x_temp;
            return;
        end
        % 选要移出基集的索引s
        [~,s_idx] = max(abs_y_B);
        s = set_B(s_idx);
        y_s = y_vect(s);
        % 计算步长向量
        temp = A(set_B,:) \ eye(n); % 等价于求逆，用左除更稳定
        t_vect = zeros(m, 1);
        t_vect(set_B_Comp) = -sign(y_s) * y_vect(set_B_Comp) .* (A(set_B_Comp, :) * temp(:, s_idx));
        % 选要加入基集的索引r
        cur_min = inf;
        cur_r = [];
        for j = set_B_Comp
            t_j = t_vect(j);
            if t_j > 1e-6 % 仅考虑正步长
                temp1 = abs(h(j)) / t_j;
                if temp1 < cur_min
                    cur_min = temp1;
                    cur_r = j;
                end
            end
        end
        if isempty(cur_r)
            error('无有效步长，问题无界');
        end
        % 基交换
        set_B = union(setdiff(set_B, s), cur_r);
        set_B_Comp = setdiff(1:m, set_B);
    end
    warning('达到最大迭代次数，返回当前解');
    x_opt = x_temp;
end

内容的提问来源于stack exchange，提问作者Shrine 420