function [cost,grad] = ricaCost(theta, x, params)

% unpack weight matrix
W = reshape(theta, params.numFeatures, params.n);

% project weights to norm ball (prevents degenerate bases)
Wold = W;
W = l2rowscaled(W, 1);

% cost function
Wx_norm = sqrt((W * x).^2 + params.epsilon);
cost = params.lambda * sum(sum(Wx_norm)) + 1/2 * sum(sum((W' * W * x - x).^2));

% gradient
% reconstruction term
grad1 = 2 * W * (W' * W * x - x) * x' + 2 * W * x * (W' * W * x - x)';

% sparsity term
grad2 = (W * x) ./ Wx_norm * x';
Wgrad = 1/2 * grad1 + params.lambda * grad2;

% unproject gradient for minFunc
grad = l2rowscaledg(Wold, W, Wgrad, 1);
grad = grad(:);

end