How to manually code 1-norm regression as a matlab function, using the below algorithm

2021-12-13

An answer to this question on Stack Overflow.

Question

I am not sure if what I have done so far is correct, and I need help with the iterative step as I don't understand what is going on in the algorithm. Here is my code. Help to finish this would be much appreciated. Thank you

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

Answer

I don't understand what's going on in your algorithm either. It's written without comments or explanations.

However, you can model your problem as a convex optimization problem. A formulation in Python using cvxpy is then quite simple and readable:

#!/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)

and gives

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