text:
function [ status, x ] = newton ( n, x0, p, f, fp, tol )

%*****************************************************************************80
%
%% NEWTON carries out the Newton iteration for the continuation system.
%
%  Discussion:
%
%    We assume that X0 is a point on the curve of solutions.
%
%    Given a starting point X, this function applies Newton's method
%    to the system of N equations:
%
%      G(X) = ( F(X)         )    <-- N-1 equations in N unknowns
%             ( X(P) - X0(P) )    <-- Hold component P fixed.
%
%    whose NxN Jacobian matrix is
%
%      GP(X) = ( FP(X)      )     <-- N-1 by N matrix.
%              ( Delta(P,J) )     <--   1 by N matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of unknowns.
%    (This should always be 2 for this program!)
%
%    Input, real X0(n), a starting point.
%
%    Input, integer P, the index of the current continuation parameter.
%
%    Input, function fx = F ( n, x ), a pointer to the
%    function that evaluates the N-1 nonlinear equations.
%
%    Input, function fpx = FP ( n, x ), a pointer to the
%    function that evaluates the N-1 by N Jacobian matrix
%    associated with the nonlinear equations.
%
%    Input, real TOL, an error tolerance.
%
%    Output, integer STATUS, the return flag.
%    0, no errors, the iteration converged.
%    1, the iteration did not converge in IT_MAX steps.
%
%    Output, real X(n), the estimated solution.
%
  verbose = 0;

  alpha = x0(p);
  x = x0;

  it = 0;
  it_max = 10;

  while ( 1 )

    if ( it_max < it )
      status = 1;
      return;
    end

    fx = f ( n, x );
    fx(n,1) = x(p) - alpha;

    fx_norm = max ( abs ( fx ) );

    if ( verbose )
      fprintf ( 1, '  %d  %g\n', it, fx_norm );
    end

    if ( fx_norm <= tol )
      status = 0;
      return;
    end

    it = it + 1;

    fpx = fp ( n, x );
    fpx(n,1:n) = 0.0;
    fpx(n,p) = 1.0;

    dx = - fpx \ fx;
    x0 = x;
    x = x0 + dx;

  end

  return
end


text:
function [nam] = get_name(dat)

%   get_name(data)  returns the name of the object
  nam=dat.algorithm.name;

%%%% SGE support
if isdeferred(dat), nam = get_name(dat.algorithm.deferred, nam, 1); end
%%%%

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function ydata = tsne_p(P, labels, no_dims)
%TSNE_P Performs symmetric t-SNE on affinity matrix P
%
%   mappedX = tsne_p(P, labels, no_dims)
%
% The function performs symmetric t-SNE on pairwise similarity matrix P
% to create a low-dimensional map of no_dims dimensions (default = 2).
% The matrix P is assumed to be symmetric, sum up to 1, and have zeros
% on the diagonal.
% The labels of the data are not used by t-SNE itself, however, they
% are used to color intermediate plots. Please provide an empty labels
% matrix [] if you don't want to plot results during the optimization.
% The low-dimensional data representation is returned in mappedX.
%
%
% (C) Laurens van der Maaten, 2010
% University of California, San Diego

    if ~exist('labels', 'var')
        labels = [];
    end
    if ~exist('no_dims', 'var') || isempty(no_dims)
        no_dims = 2;
    end

    % First check whether we already have an initial solution
    if numel(no_dims) > 1
        initial_solution = true;
        ydata = no_dims;
        no_dims = size(ydata, 2);
    else
        initial_solution = false;
    end

    % Initialize some variables
    n = size(P, 1);                                     % number of instances
    momentum = 0.5;                                     % initial momentum
    final_momentum = 0.8;                               % value to which momentum is changed
    mom_switch_iter = 250;                              % iteration at which momentum is changed
    stop_lying_iter = 100;                              % iteration at which lying about P-values is stopped
    max_iter = 1000;                                    % maximum number of iterations
    epsilon = 500;                                      % initial learning rate
    min_gain = .01;                                     % minimum gain for delta-bar-delta

    % Make sure P-vals are set properly
    P(1:n + 1:end) = 0;                                 % set diagonal to zero
    P = 0.5 * (P + P');                                 % symmetrize P-values
    P = max(P ./ sum(P(:)), realmin);                   % make sure P-values sum to one
    const = sum(P(:) .* log(P(:)));                     % constant in KL divergence
    if ~initial_solution
        P = P * 4;                                      % lie about the P-vals to find better local minima
    end

    % Initialize the solution
    if ~initial_solution
        ydata = .0001 * randn(n, no_dims);
    end
    y_incs  = zeros(size(ydata));
    gains = ones(size(ydata));

    % Run the iterations
    for iter=1:max_iter

        % Compute joint probability that point i and j are neighbors
        sum_ydata = sum(ydata .^ 2, 2);
        num = 1 ./ (1 + bsxfun(@plus, sum_ydata, bsxfun(@plus, sum_ydata', -2 * (ydata * ydata')))); % Student-t distribution
        num(1:n+1:end) = 0;                                                 % set diagonal to zero
        Q = max(num ./ sum(num(:)), realmin);                               % normalize to get probabilities

        % Compute the gradients (faster implementation)
        L = (P - Q) .* num;
        y_grads = 4 * (diag(sum(L, 1)) - L) * ydata;

        % Update the solution
        gains = (gains + .2) .* (sign(y_grads) ~= sign(y_incs)) ...         % note that the y_grads are actually -y_grads
              + (gains * .8) .* (sign(y_grads) == sign(y_incs));
        gains(gains < min_gain) = min_gain;
        y_incs = momentum * y_incs - epsilon * (gains .* y_grads);
        ydata = ydata + y_incs;
        ydata = bsxfun(@minus, ydata, mean(ydata, 1));

        % Update the momentum if necessary
        if iter == mom_switch_iter
            momentum = final_momentum;
        end
        if iter == stop_lying_iter && ~initial_solution
            P = P ./ 4;
        end

        % Print out progress
        if ~rem(iter, 10)
            cost = const - sum(P(:) .* log(Q(:)));
            disp(['Iteration ' num2str(iter) ': error is ' num2str(cost)]);
        end

        % Display scatter plot (maximally first three dimensions)
        if ~rem(iter, 10) && ~isempty(labels)
            if no_dims == 1
                scatter(ydata, ydata, 9, labels, 'filled');
            elseif no_dims == 2
                scatter(ydata(:,1), ydata(:,2), 9, labels, 'filled');
            else
                scatter3(ydata(:,1), ydata(:,2), ydata(:,3), 40, labels, 'filled');
            end
            axis tight
            axis off
            drawnow
        end
    end

repo: mobile-id


text:
function [ vX, vV, mX ] = ProjectedGdMomentum( vX, vV, hG, hP, numIterations, stepSizeGd, stepSizeMomentum )
% ----------------------------------------------------------------------------------------------- %
% [ vX, vV, mX ] = ProjectedGdMomentum( vX, vV, hG, hP, numIterations, stepSizeGd, stepSizeMomentum )
% Solves an objective function by the Projected Gradient Descent /
% Projected Gradient Method (PGM) with momentum.
% Input:
%   - vX                -   Input Vector.
%                           The starting point for the gradient iterations.
%                           Structure: Vector (numElements x 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vV                -   Momentum Vector.
%                           The buffer to hold the momentum direction
%                           vector.
%                           Structure: Vector (numElements x 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - hG                -   Gradient Function Handler.
%                           A function handler which accepts a vector of
%                           size (numElements x 1) and returns the gradient
%                           in the form of a vector (numElements x 1).
%                           Structure: Function Handler.
%                           Type: NA.
%                           Range: NA.
%   - hP                -   Gradient Function Handler.
%                           A function handler which accepts a vector of
%                           size (numElements x 1) and returns the gradient
%                           in the form of a vector (numElements x 1).
%                           Structure: Function Handler.
%                           Type: NA.
%                           Range: NA.
%   - numIterations     -   Number of Iterations.
%                           Sets the number of iterations of the gradient
%                           descent.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, ...}.
%   - stepSizeGd        -   Step Size.
%                           Sets the step size of the Gradient step.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: (0, inf).
%   - stepSizeMomentum  -   Momentum Step Size.
%                           The step size of the momentum step.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: [0, 1).
% Output:
%   - vX                -   Output Vector.
%                           The end point for the gradient iterations.
%                           Structure: Vector (numElements x 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vV                -   Acceleration Vector.
%                           The buffer to hold the momentum direction
%                           vector of last iteration. May used by a wrapper
%                           for warm start.
%                           Structure: Vector (numElements x 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - mX                -   The Path Matrix.
%                           The points of the gradient iterations.
%                           Structure: Matrix (numElements x numIterations).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References:
%   1.  Why Momentum Really Works (https://distill.pub/2017/momentum/).
% Remarks:
%   1.  B
% TODO:
%   1.  C
% Release Notes:
%   -   1.0.000     26/12/2020  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

%{
V0 = 0
V1 = beta * v0 + alpha * (1 - beta) * vG1
V2 = beta * beta * v0 + beta * alpha * (1 - beta) * vG1 + alpha * (1 -
beta) * vG2 = alpha * beta * (1 - beta) * vG1 + alpha * (1 - beta) * vG2 =
alpha * (beta * (1 - beta) * vG1 + (1 - beta) * vG2)

Then:

vW = vW - Vt = vW - alpha * Vt`

Where Vt = beta * Vt-1 + (1 - beta) * vGt

%}
stepSize    = stepSizeGd * (1 - stepSizeMomentum);
vG          = zeros(size(vX, 1), 1);

mX(:, 1) = vX;

for ii = 2:numIterations
    vG(:) = hG(vX); %<! The gradient

    % My variant
%     vV(:) = (stepSizeMoment * vV) - (stepSizeGd * vG); %<! Momentum step
%     vX(:) = hP(vX + vV); %<! Projection step

    % https://towardsdatascience.com/a84097641a5d
    % Seems to have very insignificant improvement over SGD
%     vV(:) = (stepSizeMoment * vV) + ((1 - stepSizeMoment) * vG); %<! Momentum step
%     vX(:) = hP(vX - (stepSizeGd *  vV)); %<! Projection step

    % https://ruder.io/optimizing-gradient-descent/
    vV(:) = (stepSizeMomentum * vV) + (stepSizeGd * vG); %<! Momentum step
    vX(:) = hP(vX - vV); %<! Projection step

    mX(:, ii) = vX;
end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function a = r8mat_diag_set_vector ( n, a, v )

%*****************************************************************************80
%
%% R8MAT_DIAG_SET_VECTOR sets the diagonal of a matrix to a vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of rows and columns of the matrix.
%
%    Input, real A(N,N), the matrix.
%
%    Input, real V(N), the vector to be assigned to the
%    diagonal of A.
%
%    Output, real A(N,N), the modified matrix.
%
  for i = 1 : n
    a(i,i) = v(i);
  end

  return
end

author: johannesgerer

repo: jburkardt-m