Generating function handles with large numbers of variable coefficients

4 views (last 30 days)
I've come across a number of types of optimization process which require a parameter search in a high-dimensional space where the number of coefficients to be fit is based on a dataset.
Example: performing PCA on a set of 4k images. The SVD method does not work because the matrix would be ~8 petabytes.
In this context one normalizes each picture, subtracts their average, and renormalizes these difference images. Those images are eigenvectors. One then maximizes the function
sum(a(n)*I(n)).^2
such that
sum(a(n).^2)=1
This requires a function handle of the type
y = @(b,I) b(1)*I(:,:,1)+b(2)*I(:,:,2)+ .. +b(n)*I(:,:,n);
Along with the function to be minimized,
OLS = @(b) (2-sum(sum((y(b)))).^2); % Minimize this to maximize norm(y(b))
Currently i need to write 'y' out explicitly.
Now the question: is there a way to express this in terms of matrix multiplication without explicitly writing it out?
Thank you,
Matt
  4 Comments
Matt J
Matt J on 23 May 2019
Edited: Matt J on 23 May 2019
OK, but what about my question? Have a correctly paraphrased what you are trying to do? Or if not, what is the general statement of your problem? And what are the typical sizes of I and n.
Matthew Reed
Matthew Reed on 23 May 2019
Edited: Matthew Reed on 23 May 2019
I am not currently trying to perform PCA on anything.
My problem is I had to type this in manually at one point (when I was doing PCA):
x = fmincon(@(x)(1-sum(((x(1)*probesettemp(:,1)+x(2)*probesettemp(:,2))+x(3)*probesettemp(:,3)+x(4)*probesettemp(:,4)+x(5)*probesettemp(:,5)+x(6)*probesettemp(:,6)+x(7)*probesettemp(:,7)+x(8)*probesettemp(:,8)+x(9)*probesettemp(:,9)+x(10)*probesettemp(:,10)+x(11)*probesettemp(:,11)+x(12)*probesettemp(:,12)+x(13)*probesettemp(:,13)+x(14)*probesettemp(:,14)+x(15)*probesettemp(:,15)+x(16)*probesettemp(:,16)+x(17)*probesettemp(:,17)+x(18)*probesettemp(:,18)+x(19)*probesettemp(:,19)+x(20)*probesettemp(:,20)+x(21)*probesettemp(:,21)+x(22)*probesettemp(:,22)+x(23)*probesettemp(:,23)+x(24)*probesettemp(:,24)+x(25)*probesettemp(:,25)+x(26)*probesettemp(:,26)+x(27)*probesettemp(:,27)+x(28)*probesettemp(:,28)+x(29)*probesettemp(:,29)+x(30)*probesettemp(:,30)+x(31)*probesettemp(:,31)+x(32)*probesettemp(:,32)+x(33)*probesettemp(:,33)+x(34)*probesettemp(:,34)+x(35)*probesettemp(:,35)+x(36)*probesettemp(:,36)+x(37)*probesettemp(:,37)+x(38)*probesettemp(:,38)+x(39)*probesettemp(:,39)+x(40)*probesettemp(:,40)+x(41)*probesettemp(:,41)+x(42)*probesettemp(:,42)+x(43)*probesettemp(:,43)+x(44)*probesettemp(:,44)+x(45)*probesettemp(:,45)+x(46)*probesettemp(:,46)+x(47)*probesettemp(:,47)+x(48)*probesettemp(:,48)+x(49)*probesettemp(:,49)+x(50)*probesettemp(:,50)+x(51)*probesettemp(:,51)+x(52)*probesettemp(:,52)+x(53)*probesettemp(:,53)+x(54)*probesettemp(:,54)+x(55)*probesettemp(:,55)+x(56)*probesettemp(:,56)+x(57)*probesettemp(:,57)+x(58)*probesettemp(:,58)+x(59)*probesettemp(:,59)+x(60)*probesettemp(:,60)+x(61)*probesettemp(:,61)+x(62)*probesettemp(:,62)+x(63)*probesettemp(:,63)+x(64)*probesettemp(:,64)+x(65)*probesettemp(:,65)+x(66)*probesettemp(:,66)+x(67)*probesettemp(:,67)+x(68)*probesettemp(:,68)+x(69)*probesettemp(:,69)+x(70)*probesettemp(:,70)+x(71)*probesettemp(:,71)+x(72)*probesettemp(:,72)+x(73)*probesettemp(:,73)+x(74)*probesettemp(:,74)+x(75)*probesettemp(:,75)+x(76)*probesettemp(:,76)+x(77)*probesettemp(:,77)+x(78)*probesettemp(:,78)+x(79)*probesettemp(:,79)+x(80)*probesettemp(:,80)+x(81)*probesettemp(:,81)+x(82)*probesettemp(:,82)+x(83)*probesettemp(:,83)+x(84)*probesettemp(:,84)+x(85)*probesettemp(:,85)+x(86)*probesettemp(:,86)+x(87)*probesettemp(:,87)+x(88)*probesettemp(:,88)+x(89)*probesettemp(:,89)+x(90)*probesettemp(:,90)+x(91)*probesettemp(:,91)+x(92)*probesettemp(:,92)+x(93)*probesettemp(:,93)+x(94)*probesettemp(:,94)+x(95)*probesettemp(:,95)+x(96)*probesettemp(:,96)+x(97)*probesettemp(:,97)+x(98)*probesettemp(:,98)+x(99)*probesettemp(:,99)+x(100)*probesettemp(:,100)+x(101)*probesettemp(:,101)).^2)),x0,A,b,Aeq,beq,LB,UB,nonlcon,options);
I'd like to type something like this:
x = fmincon(@(x)(1-x.*probesettemp)
but allow the contents of the array x to be a set of variational parameters, as in the equation as expressed above.

Sign in to comment.

Accepted Answer

Steven Lord
Steven Lord on 23 May 2019
You're using release R2018a, so you can take advantage of implicit expansion.
I = reshape(1:60, [3 4 5]);
b = [1 2 3 4 5];
R = reshape(b, [1, 1, numel(b)]);
S = sum(I.*R, 3);
The size of I is [3 4 5] and the size of R is [1 1 5] so the product I.*R has size [3 4 5]. Summing that product in the third dimension results in a matrix of size [3 4]. You can compare this with the result of explicitly writing out the multiplication.
S2 = b(1)*I(:,:,1)+b(2)*I(:,:,2)+ b(3)*I(:, :, 3)+b(4)*I(:, :, 4) +b(5)*I(:,:,5);
S-S2
Unfortunately the ability of the sum function to operate on multiple dimensions at a time was introduced in release R2018b, so you can't use that to simplify your OLS function. But if you could upgrade:
OLS = 2-sum(S, 'all').^2
Or combining the expressions together:
OLS = 2-sum(I.*R, 'all').^2
You don't have to define a separate variable with the reshaped b. I did that for clarity of the example.
But going back to your original question, the problem you're trying to solve is taking the pca of a matrix that's larger than can fit in memory, correct? Try storing your data as a tall array and calling pca on that tall array. The documentation for pca in release R2018a says it supports some of the syntaxes for pca on tall arrays.
  11 Comments
Matthew Reed
Matthew Reed on 24 May 2019
If you specify array operations with a function handle,
y=@(z)z*x
and then call it, with, say,
fminsearch(y,guess,opts)
Then it uses the dimensions of 'guess' to create a variable array.

Sign in to comment.

More Answers (1)

Matt J
Matt J on 23 May 2019
Edited: Matt J on 23 May 2019
I think this is what you want, but am still not totally sure. The way to calculate a sum of scalar/matrix products without explicitly writing out all n terms (if that's the question) would be,
function y=linearOp(b,I)
[p,q,n]=size(I);
y=reshape(I,[],n)*b(:);
y=reshape(y,p,q);
end

Categories

Find more on Matrices and Arrays in Help Center and File Exchange

Products


Release

R2018a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!