如何使用以下算法手动将 1 范数回归编码为 matlab 函数

Posted 2023-02-19

【中文标题】如何使用以下算法手动将 1 范数回归编码为 matlab 函数

【英文标题】：How to manually code 1-norm regression as a matlab function, using the below algorithm

【发布时间】：2022-01-07 16:08:38

【问题描述】：

我不确定到目前为止我所做的是否正确，并且我需要迭代步骤的帮助，因为我不明白算法中发生了什么。这是我的代码。帮助完成这将不胜感激。谢谢

function x_opt = out(A,b)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
b_vect = b';
m = size(A,1);
n = size(1,A);
set_B = 1:n;
set_B_Comp = n+1:m;
M = inv(A(set_B, :));
is_opt = 0;
x_temp = M*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) = -inv(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)
   x_opt = x_temp;
   ...
else
   all_index_y_vect_more_than_1 = find(abs(y_vect) >= 1);
   set_B_index_y_vect_more_than_1 = intersect(set_B, all_index_y_vect_more_than_1);
   s = set_B_index_y_vect_more_than_1(1);
   y_s = y(s)
   t_vect = zeros(m, 1);
   temp = inv(A(set_B,:));
   t_vect(set_B_Comp) = -(sign(y_s))*(y(set_B_Comp)).*(A(set_B_Comp, :)*temp(:, s));
  
   cur_min = h(set_B_Comp(1))/t_vect(set_B_Comp(1)) + 1;
   cur_r = set_B_Comp(1);
   for j = set_B_Comp
       h_j = h(j);
       t_j = t_vect(j);
       temp1 = abs(h_j)/t_j;
       if (temp1 < cur_min) && (temp1 > 0)
           cur_min = temp1;
           cur_r = j;
       end
   end
   r = cur_r;
   set_B_new = union(setdiff(set_B, s), r);
   set_B_Comp_new = setdiff(1:m,set_B_new);
   x_new = inv(A(set_B_new, :))*b_vect(set_B_new);
end
x_opt = x_temp;
end

【参考方案1】：

我也不明白你的算法发生了什么。它是在没有 cmets 或解释的情况下编写的。

但是，您可以将问题建模为凸优化问题。在 Python 中使用 cvxpy 的公式非常简单易读：

#!/usr/bin/env python3

import cvxpy as cp
import numpy as np

# Coefficient of regularization
alpha = 1e-4

# Dimensionality
N = 400
d = 20

# Synthetic data
A = np.random.randn(N, d)
b = np.random.randn(N)

# Define and solve the CVXPY problem.
x = cp.Variable(d)
objective = cp.sum_squares(A @ x - b) + alpha * cp.norm1(x)
prob = cp.Problem(cp.Minimize(objective))
optval = prob.solve()

# Print result.
print("Optimal value ", optval)
print("The optimal x is")
print(x.value)
print("The norm of the residual is ", cp.norm(A @ x - b, p=2).value)

并给予

Optimal value  328.41957961297607
The optimal x is
[-0.02041302 -0.16156503  0.07215877  0.00505087  0.01188415 -0.01029848
 -0.0237066   0.0370556   0.02205413  0.00137185  0.04055319 -0.01157271
  0.00369032  0.06349145  0.07494259 -0.04172275  0.04376864  0.02698337
 -0.04696984  0.05245699]
The norm of the residual is  18.122348149231115

【讨论】：

所以我知道库可以用来自动完成。这里的目标是手动实现这个算法。这是我正在遵循的过程：[链接] (drive.google.com/file/d/14xaCK-1Mpd8FXM-19pfFC1UTk2V9oXkQ/…)

你知道如何进行吗？