Generating function handles with large numbers of variable coefficients

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

[Matt J]

So you are trying to find the maximum singular value/vector of the matrix

A=reshape(I,[],n);

but A cannot be stored in RAM? How big is n and what is size(I)?

[Matthew Reed]

PCA is just one example of a set of optimization problems which require on-the-fly fits outputting coefficients of eigenvectors that don't lend themselves to projections. Because they're [nonlinear operators] in banach spaces.

[Matt J]

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]

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.

Answer 1

[Steven Lord]

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.

[Matthew Reed]

Thank you, Steven. That is certainly a more efficient way to compute those products than for-loops, which is my old go-to.

Would I be able to use the parameters in R as fitting parameters? Or does the fact that R has declared values keep this from working?

OLS=@(R)2-sum(I.*R, 'all').^2
params = fminsearch(OLS, guess, opts);

where params would contain the fitted values of R.

[Catalytic]

You're using release R2018a, so you can take advantage of implicit expansion.

The problem I see with that is that I.*R will consume as much memory as I itself, which is ostensibly a large array that we want to avoid copying. The method shown by Matt J avoids that, by posing the sum as an mtimes operation.

[Matthew Reed]

I seem to have confused everyone by providing an example of a context in which I'd find a compact way to generate function handles useful.

My question isn't really about handling large data sets, or manipulating arrays with assigned values. It's about writing function handles where there are a large number of variables.

[Catalytic]

My question isn't really about handling large data sets, or manipulating arrays with assigned values. It's about writing function handles where there are a large number of variables.

Matt and Steve's responses do answer that question, though.The answer is not to write such a function at all, but rather to vectorize your operations, writing them in terms of a single vector variable instead of lots of separate scalar variables.

[Matthew Reed]

I must be getting confused by the lack of @ in any of the code that I'm reading.

If i want to solve an optimization problem, and I write down the routine as:

y=@(a)a.*X
params = fminsearch(OLS, guess, opts);

it should just work? Does the size of 'guess' define the number of elements of 'a'? Or does 'a' have to be in the workspace, containing any set of numbers, as long as it's shaped correctly?

[Catalytic]

I guess you figured it out, since you Accept-clicked an answer?

[madhan ravi]

@Catalytic: This is not my business but just curious. Are you Matt's workmate by any chance?

[Catalytic]

@madhan ravi: Matt Reed or Matt J? I am the workmate of neither, but why do you ask?

[madhan ravi]

No worries, was just curious.

[Matt J]

@Catalytic: where do you work, if you don't mind my asking?

[Matthew Reed]

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.

Answer 2

[Matt J]

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