Kappa and fista (#14)

* added GD kappa effect & Image deblurring w/ FISTA

* added GD kappa effect & Image deblurring w/ FISTA

* Update README.md

* Update README.md

* Update README.md

* Rename Image_Related/Image Deblurring/generic_grad.m to Image_Related/Image_Deblurring/generic_grad.m

* Rename Image_Related/Image Deblurring/fista.m to Image_Related/Image_Deblurring/fista.m

* Rename Image_Related/Image Deblurring/exact_quad.m to Image_Related/Image_Deblurring/exact_quad.m

* Rename Image_Related/Image Deblurring/const_step.m to Image_Related/Image_Deblurring/const_step.m

* Rename Image_Related/Image Deblurring/blur.m to Image_Related/Image_Deblurring/blur.m

* Rename Image_Related/Image Deblurring/README.md to Image_Related/Image_Deblurring/README.md

* Rename Image_Related/Image Deblurring/GD_vs_FISTA_wrapper.m to Image_Related/Image_Deblurring/GD_vs_FISTA_wrapper.m

* Rename

Image_Related/Image_Deblurring/GD_vs_FISTA_wrapper.m

@@ -0,0 +1,98 @@
+% Basic Example
+% sweden = [1,2,1,1,1;1,2,1,1,1;2,2,2,2,2;1,2,1,1,1;1,2,1,1,1];
+% colormap(parula)
+% imagesc(sweden)
+
+
+% Display the original and blurred images
+% colormap ( gray )
+% imagesc ( reshape (x , 256 , 256 ) )
+% figure()
+% colormap ( gray )
+% imagesc ( reshape (b , 256 , 256 ) )
+[A, b , x ] = blur ( 256 , 5 , 1 ) ;
+
+f = @(x) (norm(A*x-b))^2;
+
+% Cmputing gradient elements beforehand
+Y1 =A'*A;
+Y = 2*Y1;
+Z = 2*A'*b;
+g = @(x) Y*x-Z; %2*A'*(A*x-b)
+
+x0 = zeros(size(x));
+% tic;
+% % ***  Part 1  ***
+% [~ ,x_1,x_10,x_100,x_1000] = generic_grad_q3(f, g, exact_quad_q3(Y), x0,1000);
+% disp('Generic Grad with Exact Quad complete. Time taken:')
+% toc
+% figure('Name','Image after 1 Iteration of GD with Exact Quad');
+% colormap ( gray )
+% imagesc ( reshape (x_1 , 256 , 256 ) )
+% 
+% figure('Name','Image after 10 Iterations of GD with Exact Quad');
+% colormap ( gray )
+% imagesc ( reshape (x_10 , 256 , 256 ) )
+% 
+% figure('Name','Image after 100 Iterations of GD with Exact Quad');
+% colormap ( gray )
+% imagesc ( reshape (x_100 , 256 , 256 ) )
+% 
+% figure('Name','Image after 1000 Iterations of GD with Exact Quad');
+% colormap ( gray )
+% imagesc ( reshape (x_1000 , 256 , 256 ) )
+% 
+% %***  Part 2  ***
+L = 2*eigs(Y1,1);
+% tic;
+% [~ ,x_1,x_10,x_100,x_1000] = generic_grad_q3(f, g, const_step_q3(1/L), x0, 1000);
+% disp('Generic Grad with const step complete. Time taken:')
+% toc
+% figure('Name','Image after 1 Iteration of GD with Const Step');
+% colormap ( gray )
+% imagesc ( reshape (x_1 , 256 , 256 ) )
+% 
+% figure('Name','Image after 10 Iterations of GD with Const Step');
+% colormap ( gray )
+% imagesc ( reshape (x_10 , 256 , 256 ) )
+% 
+% figure('Name','Image after 100 Iterations of GD with Const Step');
+% colormap ( gray )
+% imagesc ( reshape (x_100 , 256 , 256 ) )
+% 
+% figure('Name','Image after 1000 Iterations of GD with Const Step');
+% colormap ( gray )
+% imagesc ( reshape (x_1000 , 256 , 256 ) )
+
+% ***  Part 3  ***
+% tic;
+% [~ ,x_1,x_10,x_100,x_1000] = fista(f, g, L, x0, 1000);
+% toc
+% figure('Name','Image after 1 Iteration of FISTA');
+% colormap ( gray )
+% imagesc ( reshape (x_1 , 256 , 256 ) )
+% 
+% figure('Name','Image after 10 Iteration of FISTA');
+% colormap ( gray )
+% imagesc ( reshape (x_10 , 256 , 256 ) )
+% 
+% figure('Name','Image after 100 Iteration of FISTA');
+% colormap ( gray )
+% imagesc ( reshape (x_100 , 256 , 256 ) )
+% 
+% figure('Name','Image after 1000 Iteration of FISTA');
+% colormap ( gray )
+% imagesc ( reshape (x_1000 , 256 , 256 ) )
+
+% Plot gs versus Time
+eps = 10^(-5);
+[x_const,fs_const,gs_const,ts_const] = generic_grad(f, g, const_step(1/L), x0, 1000);
+[x_fista,fs_fista,gs_fista,ts_fista] = fista(f, g, L, x0, eps);
+figure('Name', 'f vs Time Comparison');
+semilogy(ts_const, fs_const);
+hold on;
+semilogy(ts_fista, gs_fista);
+title('Semilog Plot of f vs Time');
+xlabel('Time') ;
+ylabel('f(x)') ;
+legend({'const','fista'},'Location','northeast');

Image_Related/Image_Deblurring/README.md

@@ -0,0 +1,3 @@
+[Results](/Image_Related/Image_Deblurring/doc/GD_vs_FISTA_wrapper.pdf)
+
+To reconstruct run [GD_vs_FISTA_wrapper.m](/Image_Related/Image_Deblurring/GD_vs_FISTA_wrapper.m)

Image_Related/Image_Deblurring/blur.m

@@ -0,0 +1,83 @@
+function [A,b,x] = blur(N,band,sigma) 
+%BLUR Test problem: digital image deblurring. 
+% 
+% function [A,b,x] = blur(N,band,sigma) 
+% 
+% The matrix A is an N*N-by-N*N symmetric, doubly block Toeplitz matrix that 
+% models blurring of an N-by-N image by a Gaussian point spread function. 
+% It is stored in sparse matrix format. 
+% 
+% In each Toeplitz block, only matrix elements within a distance band-1 
+% from the diagonal are nonzero (i.e., band is the half-bandwidth). 
+% If band is not specified, band = 3 is used. 
+% 
+% The parameter sigma controls the width of the Gaussian point spread 
+% function and thus the amount of smoothing (the larger the sigma, the wider 
+% the function and the more ill posed the problem).  If sigma is not 
+% specified, sigma = 0.7 is used. 
+% 
+% The vector x is a columnwise stacked version of a simple test image, while 
+% b holds a columnwise stacked version of the blurrred image; i.e, b = A*x. 
+
+% Per Christian Hansen, IMM, 11/11/97. 
+
+% Initialization. 
+if (nargin < 2), band = 3; end 
+band = min(band,N); 
+if (nargin < 3), sigma = 0.7; end 
+
+% Construct the matrix as a Kronecker product. 
+z = [exp(-([0:band-1].^2)/(2*sigma^2)),zeros(1,N-band)]; 
+A = toeplitz(z); 
+A = sparse(A); 
+A = (1/(2*pi*sigma^2))*kron(A,A); 
+
+% Generate x and b, if required. 
+if (nargout > 1) 
+
+  % Start with an image of all zeros. 
+  x = zeros(N,N); 
+  N2 = round(N/2); 
+  N3= round(N/3); 
+  N6 = round(N/6); 
+  N12 = round(N/12); 
+
+  % Add a large ellipse. 
+  T = zeros(N6,N3); 
+  for i=1:N6, for j=1:N3 
+    if ( (i/N6)^2 + (j/N3)^2 < 1 ), T(i,j) = 1; end 
+  end, end 
+  T = [fliplr(T),T]; 
+  T = [flipud(T);T]; 
+  x(2+[1:2*N6],N3-1+[1:2*N3]) = T; 
+
+  % Add a smaller ellipse. 
+  T = zeros(N6,N3); 
+  for i=1:N6, for j=1:N3 
+    if ( (i/N6)^2 + (j/N3)^2 < 0.6 ), T(i,j) = 1; end 
+  end, end 
+  T = [fliplr(T),T]; 
+  T = [flipud(T);T]; 
+  x(N6+[1:2*N6],N3-1+[1:2*N3]) = x(N6+[1:2*N6],N3-1+[1:2*N3]) + 2*T; 
+  % Correct for overlap. 
+  f = find(x==3); 
+  x(f) = 2*ones(size(f)); 
+
+  % Add a triangle. 
+  T = triu(ones(N3,N3)); 
+  [mT,nT] = size(T); 
+  x(N3+N12+[1:nT],1+[1:mT]) = 3*T; 
+
+  % Add a cross. 
+  T = zeros(2*N6+1,2*N6+1); 
+  [mT,nT] = size(T); 
+  T(N6+1,1:nT) = ones(1,nT); 
+  T(1:mT,N6+1) = ones(mT,1); 
+  x(N2+N12+[1:mT],N2+[1:nT]) = 4*T; 
+
+  % Make sure x is N-times-N, and stack the columns of x. 
+  x = reshape(x(1:N,1:N),N^2,1); 
+
+  % Compute the blurred image. 
+  b = A*x; 
+end 

Image_Related/Image_Deblurring/const_step.m

@@ -0,0 +1,10 @@
+function lsearch = const_step(s)
+
+% verify s is a real positive scalar
+if s<=0 || ~isreal(s) || ~isscalar(s)
+        error("s must be a real positive scalar")
+end
+
+lsearch =@(f,x_k,g_k) s;
+end
+

Image_Related/Image_Deblurring/exact_quad.m

@@ -0,0 +1,4 @@
+function lsearch = exact_quad_q3(A)
+lsearch = @(f,xk,gk) ( gk'*gk )/(gk'*A*gk); % leading matrix is (1/2)*A
+end
+

Image_Related/Image_Deblurring/fista.m

@@ -0,0 +1,44 @@
+function [x, fs, gs, ts] = fista(f, gf, L, x0, eps)
+
+% Input checks
+if ~isa(f,'function_handle') || ~isa(gf,'function_handle') 
+	error("f and gf must be functions")
+end
+if eps<=0 || ~isreal(eps) || ~isscalar(eps)
+        error("eps must be a real positive scalar")
+end
+if ~isreal(x0)
+    error("x0 must be real")
+end
+
+% Start Timer
+tic
+
+% Initialize Paramaters
+x = x0;
+y = x0;
+grad = gf(x);
+grad_norm = norm(grad);
+
+% Result containers
+gs = grad_norm;
+fs = f(x0);
+ts = 0;
+t_k = 1;
+
+ while (grad_norm >= eps) 
+    old_x = x;
+    x = y-((1/L )*grad);
+    old_t_k = t_k;
+    t_k = (1+sqrt(1+4*((old_t_k)^2)))/2;
+    y = x+((old_t_k-1)/t_k)*(x-old_x);
+    grad = gf(y);
+    grad_norm = norm(grad);
+
+    % Recording interim results
+    fs = [fs, f(y)];
+    gs = [gs, grad_norm];
+    ts = [ts, toc]; 
+ end
+end
+

Image_Related/Image_Deblurring/generic_grad.m

@@ -0,0 +1,52 @@
+function [x ,x_1,x_10,x_100,x_1000] = generic_grad_q3(f , gf , lsearch , x0 , MAX_ITERATIONS)
+
+% Input checks
+if ~isa(f,'function_handle') || ~isa(gf,'function_handle') || ~isa(lsearch,'function_handle')
+	error("f, gf and lsearch  must be functions")
+end
+if ~isscalar(MAX_ITERATIONS) || ~(MAX_ITERATIONS==floor(MAX_ITERATIONS)) || ~isreal(MAX_ITERATIONS)
+    error("MAX_ITERATIONS must be a real scalar integer")
+end
+if ~isreal(x0)
+    error("x0 must be real")
+end
+
+% Initialize Paramaters
+x = x0;
+grad = gf(x);
+% Result containers
+
+iteration = 0;
+
+% Execute Generic Gradient Descent Algorithm
+% while (grad_norm > eps && iteration < MAX_ITERATIONS)
+while (iteration < MAX_ITERATIONS)
+
+    t_k = lsearch(f,x,grad);
+    x = x - (t_k * grad);
+    grad = gf(x);
+
+    % Recording interim results
+    iteration=iteration+1;
+    if(iteration == 1)
+        x_1=x;
+        disp(['Iteration ', num2str(iteration), ' complete'])
+    end
+    if(iteration == 10)
+        x_10=x;
+        disp(['Iteration ', num2str(iteration), ' complete'])
+    end
+    if(iteration == 100)
+     x_100=x;
+    end
+    if(iteration == 1000)
+     x_1000=x;
+    end
+
+    if(mod(iteration,100) == 0)
+        disp(['Iteration ', num2str(iteration), ' complete'])
+    end
+
+end
+end
+

Optimization/Gradient_Descent_kappa_coef_effect/Gradient_Descent_kappa_influence.m

@@ -0,0 +1,38 @@
+% Parameter Initialization
+x0 = [1000;1000;1000];
+
+F=[2,2,3; 2,5,5; 3,5,7];
+G=[10,8,9;8,13,11;9,11,15];
+
+g=@(x) 0.5*x'*G*x;
+f=@(x) 0.5*x'*F*x;
+
+grad_g=@(x) G*x;
+grad_f=@(x) F*x;
+
+step_G = 1/eigs(G,1);
+step_F = 1/eigs(F,1);
+
+K_G = cond(G)
+K_F = cond(F)
+k_ratio = K_F/K_G
+
+disp('Generic Grad for f complete. ')
+[~, fs_f, iterations_f]  = generic_grad(f, grad_f, const_step(step_F), x0, 100000, 10^-5);
+
+disp('Generic Grad for g complete.')
+[~, fs_g, iterations_g]= generic_grad(g, grad_g, const_step(step_G), x0, 100000, 10^-5);
+
+% Plot fs versus Number of Iterations
+figure('Name','fs vs Iters Comparison');
+semilogy(fs_f);
+hold on;
+semilogy(fs_g);
+
+title('Semilogy Plot of fs & gs vs  Number of Iterations');
+xlabel('Iteration');
+ylabel('function value');
+legend({'fs\_f','fs\_g'},'Location','southwest');
+
+iterations_ratio = iterations_f/iterations_g
+

Optimization/Gradient_Descent_kappa_coef_effect/README.md

@@ -0,0 +1,16 @@
+Sample Output:
+```matlab
+K_G =
+   10.9129
+K_F =
+   25.4916
+k_ratio =
+    2.3359
+Generic Grad for f complete. 
+Total iterations taken: 250
+Generic Grad for g complete.
+Total iterations taken: 105
+iterations_ratio =
+    2.3810
+```
+![graph](/Optimization/Gradient_Descent_kappa_coef_effect/sample_output.jpg)

Optimization/Gradient_Descent_kappa_coef_effect/const_step.m

@@ -0,0 +1,10 @@
+function lsearch = const_step(s)
+
+% verify s is a real positive scalar
+if s<=0 || ~isreal(s) || ~isscalar(s)
+        error("s must be a real positive scalar")
+end
+
+lsearch =@(f,x_k,g_k) s;
+end
+