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

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

Richard · Accepted Answer · 2021-12-13 02:52:16Z

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

So I know libraries can be used to do it automatically. The goal here is to manually implement this algorithm. This is the process I am following: [link] (drive.google.com/file/d/14xaCK-1Mpd8FXM-19pfFC1UTk2V9oXkQ/…)

Collectives™ on Stack Overflow

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

1 Answer 1

2 Comments

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

2 Comments

Related