%% Jacobian Multiply Function with Linear Least Squares
% Using a Jacobian multiply function, you can solve a least-squares problem of the form
% \[\mathop {\min }\limits_x \frac{1}{2}\left\| {C \cdot x - d} \right\|_2^2\]
% such that lb ≤ x ≤ ub, for problems where C is very large, perhaps too large to be stored. For this technique, use the 'trust-region-reflective' algorithm.
% For example, consider a problem where C is a 2n-by-n matrix based on a circulant matrix. The rows of C are shifts of a row vector v. This example has the row vector v with elements of the form (–1)^{k+1}/k:
% v = [1,-1/2,1/3,-1/4,\ldots,-1/n],
% where the elements are cyclically shifted.
% \[C = \left[ {\begin{array}{*{20}{c}}
% 1&{ - 1/2}&{1/3}&{...}&{ - 1/n}\\
% { - 1/n}&1&{ - 1/2}&{...}&{1/(n - 1)}\\
% {1/(n - 1)}&{ - 1/n}&1&{...}&{ - 1/(n - 2)}\\
%  \vdots & \vdots & \vdots & \ddots & \vdots \\
% { - 1/2}&{1/3}&{ - 1/4}&{...}&1\\
% 1&{ - 1/2}&{1/3}&{...}&{ - 1/n}\\
% { - 1/n}&1&{ - 1/2}&{...}&{1/(n - 1)}\\
% {1/(n - 1)}&{ - 1/n}&1&{...}&{ - 1/(n - 2)}\\
%  \vdots & \vdots & \vdots & \ddots & \vdots \\
% { - 1/2}&{1/3}&{ - 1/4}&{...}&1
% \end{array}} \right].\]
% This least-squares example considers the problem where
% d = [n-1,n-2,\ldots,-n],
% and the constraints are -5\le x_i\le 5 for i = 1,\ldots,n.
% For large enough n, the dense matrix C does not fit into computer memory (n=10,000 is too large on one tested system).
% A Jacobian multiply function has the following syntax.
% w = jmfcn(Jinfo,Y,flag)
% Jinfo is a matrix the same size as C, used as a preconditioner. If C is too large to fit into memory, Jinfo should be sparse. Y is a vector or matrix sized so that C*Y or C'*Y works as matrix multiplication. flag tells jmfcn which product to form:
% flag > 0 ⇒  w = C*Y
% flag < 0 ⇒  w = C'*Y
% flag = 0 ⇒  w = C'*C*Y
% Because C is such a simply structured matrix, you can easily write a Jacobian multiply function in terms of the vector v, without forming C. Each row of C*Y is the product of a circularly shifted version of v times Y. Use circshift to circularly shift v.
% To compute C*Y, compute v*Y to find the first row, then shift v and compute the second row, and so on.
% To compute C'*Y, perform the same computation, but use a shifted version of temp, the vector formed from the first row of C':
% temp = [fliplr(v),fliplr(v)];
% temp = [circshift(temp,1,2),circshift(temp,1,2)]; % Now temp = C'(1,:)
% To compute C'*C*Y, simply compute C*Y using shifts of v, and then compute C' times the result using shifts of fliplr(v).
% The helper function lsqcirculant3 is a Jacobian multiply function that implements this procedure; it appears at the end of this example.
% The dolsqJac3 helper function at the end of this example sets up the vector v and calls the solver lsqlin using the lsqcirculant3 Jacobian multiply function.
% When n = 3000, C is an 18,000,000-element dense matrix. Determine the results of the dolsqJac3 function for n = 3000 at selected values of x, and display the output structure.

[x,resnorm,residual,exitflag,output] = dolsqJac3(3000);
disp(x(1))
disp(x(1500))
disp(x(3000))
disp(output)

%
%% Helper Functions
% This code creates the lsqcirculant3 helper function.

function w = lsqcirculant3(Jinfo,Y,flag,v)
% This function computes the Jacobian multiply function
% for a 2n-by-n circulant matrix example.

if flag > 0
    w = Jpositive(Y);
elseif flag < 0
    w = Jnegative(Y);
else
    w = Jnegative(Jpositive(Y));
end

    function a = Jpositive(q)
        % Calculate C*q
        temp = v;

        a = zeros(size(q)); % Allocating the matrix a
        a = [a;a]; % The result is twice as tall as the input.

        for r = 1:size(a,1)
            a(r,:) = temp*q; % Compute the rth row
            temp = circshift(temp,1,2); % Shift the circulant
        end
    end

    function a = Jnegative(q)
        % Calculate C'*q
        temp = fliplr(v);
        temp = circshift(temp,1,2); % Shift the circulant for C'

        len = size(q,1)/2; % The returned vector is half as long
        % as the input vector.
        a = zeros(len,size(q,2)); % Allocating the matrix a

        for r = 1:len
            a(r,:) = [temp,temp]*q; % Compute the rth row
            temp = circshift(temp,1,2); % Shift the circulant
        end
    end
end

% This code creates the dolsqJac3 helper function.

function [x,resnorm,residual,exitflag,output] = dolsqJac3(n)
%
r = 1:n-1; % Index for making vectors

v(n) = (-1)^(n+1)/n; % Allocating the vector v
v(r) =( -1).^(r+1)./r;

% Now C should be a 2n-by-n circulant matrix based on v,
% but it might be too large to fit into memory.

r = 1:2*n;
d(r) = n-r;

Jinfo = [speye(n);speye(n)]; % Sparse matrix for preconditioning
% This matrix is a required input for the solver;
% preconditioning is not used in this example.

% Pass the vector v so that it does not need to be
% computed in the Jacobian multiply function.
options = optimoptions('lsqlin','Algorithm','trust-region-reflective',...
    'JacobianMultiplyFcn',@(Jinfo,Y,flag)lsqcirculant3(Jinfo,Y,flag,v));

lb = -5*ones(1,n);
ub = 5*ones(1,n);

[x,resnorm,residual,exitflag,output] = ...
    lsqlin(Jinfo,d,[],[],[],[],lb,ub,[],options);
end

% Copyright 2019 The MathWorks, Inc.