Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
fminsearch for functions

Subject: fminsearch for functions

From: Karl

Date: 19 Sep, 2011 20:03:27

Message: 1 of 12

Hi, I have a problem and I was wondering if I could use fminsearch (or some related function) to help me out. Essentially I have 15 unknown functions which are combined in some particular way (I know how they are combined, if I have the 15 I can easily enough combine them) to give me a measured matrix (I know this matrix). So my problem is essentially I want to minimze the norm of (the measured matrix minus the "what i predict" matrix from 15 functions). It is essentailly a minimzation of multiple functions problem, whereas a typical minimization problem says "find the minimum of a function" i'm saying "find the minimum of the norm of (1 matix - the other matrix)" by fiddling with 15 functions (vectors). Is this possible? I have considered doing simple-minded iterative techniques myself but I believe they will be too computationally intensive to run in the amount of time I desire, thus I
was hoping some MATLAB iterative technique could help me out. Thank you very much for any help.

Subject: fminsearch for functions

From: Matt J

Date: 19 Sep, 2011 20:43:26

Message: 2 of 12

"Karl " <karllandheer@gmail.com> wrote in message <j5876f$g8s$1@newscl01ah.mathworks.com>...
>
 It is essentailly a minimzation of multiple functions problem, whereas a typical minimization problem says "find the minimum of a function" i'm saying "find the minimum of the norm of (1 matix - the other matrix)" by fiddling with 15 functions
===============

This seems like a false distinction. The norm of the difference of two matrices is a single function and you're proposing to minimize that. So, the 2nd case you describe is a typical minimization problem as well, just like the first case.

You can try FMINSEARCH, but the functions in the Optimization Toolbox will probably be better, particularly considering the number of unknown variables that you have.

Subject: fminsearch for functions

From: Karl

Date: 21 Sep, 2011 16:41:27

Message: 3 of 12

You are right, I didn't realize fminsearch was able to accept matrices as an input, thanks for this. I am, however, having another problem. Here is my function and I wish to minimize the L2Norm:

function [L2Norm] = forwardScatterProjection(radProf)

load('GeneratedPEMultiBeam_noisyBeams')
numberOfRows = size(RawImage,1);
numberOfColumns = size(RawImage,2);
numberOfBeams = 5;

maxRadiusIndex = [(j0) (616-j0)];
maxRadiusIndex = min(maxRadiusIndex);

centreLocationY = 305; %all y locations should be the same
% centreLocationX = 349:35:839; %there is 15 beams so 15 X location centres
centreLocationX = positionLeft:pixelSpacing:positionRight; %there is 15 beams so 15 X location centres
%centreLocationX = 580:1:635; %there is 15 beams so 15 X location centres
reconstructedMatrix = zeros(numberOfRows,numberOfColumns);
radiusIndex = 1:maxRadiusIndex; %the size of the radii for the angles
centreMatrix = zeros(numberOfRows,numberOfColumns);


for beamNumber = 1:numberOfBeams
    centreMatrix(centreLocationY,centreLocationX(beamNumber)) = 1;
    distanceMatrixCell{beamNumber,:,:} = floor(bwdist(centreMatrix));
    centreMatrix(centreLocationY,centreLocationX(beamNumber)) = 0;
end

reconstructedMatrix = zeros(NumberOfRows,NumberOfColumns);%initialize it to zeros then rebuild it
%forward project
for beamNumber = 1:numberOfBeams
    distanceMatrix = distanceMatrixCell{beamNumber}; %matrix which shows the pixel distance to the centre
    for ii = 1:numberOfRows
        for jj = 1:numberOfColumns
            if (distanceMatrix(ii,jj) > 0 && distanceMatrix(ii,jj) < maxRadiusIndex && (ii < Rod(1) || ii > Rod(RodCounter)))
                reconstructedMatrix(ii,jj) = reconstructedMatrix(ii,jj) + radProf(beamNumber,distanceMatrix(ii,jj));
            end
        end
    end
end

figure;imagesc(reconstructedMatrix)
L2Norm = norm(reconstructedMatrix-RawImage)

I call my function in the following manner:
radProf = fminsearch('forwardScatterProjection', ones(5,305))

And the problem I seem to be having is that it does not seem to be attempting to minimize by iterating on the matrix "radProf". Each iteration it spits out the exact same L2Norm. I'm not sure what I'm doing wrong, my program for whatever reason not compatible with the minimizatino algorithm used (seems unlikely) or what? I have tried it with simpler functions in its place and it works but I'm not certain what I'm doing wrong. Thanks again for your continued help.

"Matt J" wrote in message <j589hd$pr8$1@newscl01ah.mathworks.com>...
> "Karl " <karllandheer@gmail.com> wrote in message <j5876f$g8s$1@newscl01ah.mathworks.com>...
> >
> It is essentailly a minimzation of multiple functions problem, whereas a typical minimization problem says "find the minimum of a function" i'm saying "find the minimum of the norm of (1 matix - the other matrix)" by fiddling with 15 functions
> ===============
>
> This seems like a false distinction. The norm of the difference of two matrices is a single function and you're proposing to minimize that. So, the 2nd case you describe is a typical minimization problem as well, just like the first case.
>
> You can try FMINSEARCH, but the functions in the Optimization Toolbox will probably be better, particularly considering the number of unknown variables that you have.

Subject: fminsearch for functions

From: Matt J

Date: 21 Sep, 2011 17:42:28

Message: 4 of 12

"Karl " <karllandheer@gmail.com> wrote in message <j5d43n$je4$1@newscl01ah.mathworks.com>...
>
> I call my function in the following manner:
> radProf = fminsearch('forwardScatterProjection', ones(5,305))
=============

So, you have 5*305=1525 unknown variables? FMINSEARCH is unlikely to do well with such a large number of unknowns. It can be barely handle 6 or 7.



> Each iteration it spits out the exact same L2Norm.

Are you sure that forwardScatterProjection is capable of returning more than one value? What tests have you done to make sure it is giving correct output?

Subject: fminsearch for functions

From: Karl

Date: 21 Sep, 2011 17:59:29

Message: 5 of 12

"Matt J" wrote in message <j5d7m4$3ud$1@newscl01ah.mathworks.com>...
> "Karl " <karllandheer@gmail.com> wrote in message <j5d43n$je4$1@newscl01ah.mathworks.com>...
> >
> > I call my function in the following manner:
> > radProf = fminsearch('forwardScatterProjection', ones(5,305))
> =============
>
> So, you have 5*305=1525 unknown variables? FMINSEARCH is unlikely to do well with such a large number of unknowns. It can be barely handle 6 or 7.
>
>
>
> > Each iteration it spits out the exact same L2Norm.
>
> Are you sure that forwardScatterProjection is capable of returning more than one value? What tests have you done to make sure it is giving correct output?

Yes, I do have 1525 unknown variables. This forwardScatterProjection algorithm is taken from a piece of code I wrote which solves the problem using a MLEM-based iterative method (modified version of a CT-reconstruction algorithm). However this algorithm seems to fall apart when I add in a 6th beam (and I was ideally hoping to do 15 beams). But this method did work perfectly for 5, so I know that the forwardProjectionScatter is doing what I hope it to do (i've also tested it by looking at the output of reconstructedMatrix from various starting points, and it output what I was expecting. e.g., if i give it ones(5,305) it will give me five circles of value 1 centred at different locations and superimposed).

Is there any tool in sftool (I don't have it on my version of matlab at work but I could upgrade if need be) that would be able to handle this large number of unknowns or should I be pleased that I even have an algorithm that can handle 5 beams?

Subject: fminsearch for functions

From: Matt J

Date: 21 Sep, 2011 18:20:26

Message: 6 of 12

"Karl " <karllandheer@gmail.com> wrote in message <j5d8m1$7a6$1@newscl01ah.mathworks.com>...
>
> Is there any tool in sftool (I don't have it on my version of matlab at work but I could upgrade if need be) that would be able to handle this large number of unknowns or should I be pleased that I even have an algorithm that can handle 5 beams?
=============

I wouldn't know the relevance to this of sftool. Like I told you before, your best bet would probably be the Optimization Toolbox. On the other hand, I don't see why MLEM would fail, except perhaps if by adding the additional beam, you've somehow made the problem ill-posed.

If you describe the problem mathematically, rather than in code, maybe it will be clear enough to the community to propose better customized algorithms.

Subject: fminsearch for functions

From: Karl

Date: 21 Sep, 2011 18:56:29

Message: 7 of 12

"Matt J" wrote in message <j5d9ta$che$1@newscl01ah.mathworks.com>...
> "Karl " <karllandheer@gmail.com> wrote in message <j5d8m1$7a6$1@newscl01ah.mathworks.com>...
> >
> > Is there any tool in sftool (I don't have it on my version of matlab at work but I could upgrade if need be) that would be able to handle this large number of unknowns or should I be pleased that I even have an algorithm that can handle 5 beams?
> =============
>
> I wouldn't know the relevance to this of sftool. Like I told you before, your best bet would probably be the Optimization Toolbox. On the other hand, I don't see why MLEM would fail, except perhaps if by adding the additional beam, you've somehow made the problem ill-posed.
>
> If you describe the problem mathematically, rather than in code, maybe it will be clear enough to the community to propose better customized algorithms.

OK my problem is simple enough. I have a matrix, call it M which is essentially generated by the superposition of some number of (call it n) "radial profiles" called radProf. The way this superposition takes place is given by the code above, but it's essentailly taking the values from the function and making them circulary distributed. So for example radProf(1) would be the centre, and radProf(2) would be the value one pixel away from the centre in any direction, etc. My problem is to recover these individual radial profiles frmo the measured Matrix. Currently this "measured" matrix is being generated by a program where I give it the radial profiles and then I run my code to see "how well" it recovers the given radial profiles. To write it down mathematically the problem is as follows

A[p1, p2 ... p15] = p1*s1 + p2*s2 + .... p15*s15 = M

where s1, s2, ... s15 are "shifting matrices" which essentially take into accuont the fact that all the radial profiles create circles with different centres. And p1, p2, and p15 are the "forward projection" of the radial profiles (i.e., take the function and generate the circlular pattern). A is simply some operator which is defined by the first equal sign.

I'm not sure why the MLEM code breaks down, I can post it if someone would like to see if there's an error but the fact that it works perfectly for 1-5 beams and then for 6 it gives a strange answer (close to what I input but not quite). Anyways I hope this is clear but it's difficult to explain here. Thank you very much for your help I really appreciate it.

Subject: fminsearch for functions

From: Matt J

Date: 21 Sep, 2011 19:40:33

Message: 8 of 12

"Karl " <karllandheer@gmail.com> wrote in message <j5dc0t$ljn$1@newscl01ah.mathworks.com>...
>
>To write it down mathematically the problem is as follows
>
> A[p1, p2 ... p15] = p1*s1 + p2*s2 + .... p15*s15 = M
>
> where s1, s2, ... s15 are "shifting matrices" which essentially take into accuont the fact that all the radial profiles create circles with different centres. And p1, p2, and p15 are the "forward projection" of the radial profiles (i.e., take the function and generate the circlular pattern). A is simply some operator which is defined by the first equal sign.
===================


First of all, are you sure you want this as your comparison function

 L2Norm = norm(reconstructedMatrix-RawImage)

or do you really want this
 
  L2Norm = norm(reconstructedMatrix(:)-RawImage(:))

Can you see why they are different? If you want the latter, then your minimization problem is equivalent to solving some linear system

  Q*radProf(:)=M(:)

for some matrix Q and you can just solve it via

 radProf(:)=Q\M(:)

If Q is too big to generate in explicit form (even when sparse), then you can use Matlab's PCG which will let you specify the operation Q*radProf(:) in functional form.



> I'm not sure why the MLEM code breaks down,

I don't see why you would be applying MLEM if you're doing L2 norm minimization. You said you were adapting MLEM from emission tomography (you said CT but I think you really meant emission tomography). Either way, MLEM applies to Poisson loglikelihood minimization, not least squares.

Subject: fminsearch for functions

From: Karl

Date: 21 Sep, 2011 20:08:26

Message: 9 of 12

"Matt J" wrote in message <j5dejh$2nr$1@newscl01ah.mathworks.com>...
> "Karl " <karllandheer@gmail.com> wrote in message <j5dc0t$ljn$1@newscl01ah.mathworks.com>...
> >
> >To write it down mathematically the problem is as follows
> >
> > A[p1, p2 ... p15] = p1*s1 + p2*s2 + .... p15*s15 = M
> >
> > where s1, s2, ... s15 are "shifting matrices" which essentially take into accuont the fact that all the radial profiles create circles with different centres. And p1, p2, and p15 are the "forward projection" of the radial profiles (i.e., take the function and generate the circlular pattern). A is simply some operator which is defined by the first equal sign.
> ===================
>
>
> First of all, are you sure you want this as your comparison function
>
> L2Norm = norm(reconstructedMatrix-RawImage)
>
> or do you really want this
>
> L2Norm = norm(reconstructedMatrix(:)-RawImage(:))
>
> Can you see why they are different? If you want the latter, then your minimization problem is equivalent to solving some linear system
>
> Q*radProf(:)=M(:)
>
> for some matrix Q and you can just solve it via
>
> radProf(:)=Q\M(:)
>
> If Q is too big to generate in explicit form (even when sparse), then you can use Matlab's PCG which will let you specify the operation Q*radProf(:) in functional form.
>
>
>
> > I'm not sure why the MLEM code breaks down,
>
> I don't see why you would be applying MLEM if you're doing L2 norm minimization. You said you were adapting MLEM from emission tomography (you said CT but I think you really meant emission tomography). Either way, MLEM applies to Poisson loglikelihood minimization, not least squares.

OK in my MLEM code I am just using the L2Norm as a kind of metric to show how close my reconstructed matrix is to my original matrix (I want it to be 0). I realize that my code will not, in general, find the minimal value of L2Norm but, correct me if i'm wrong, the most likely solution will be the one that minimizes the least squares IF there is a perfect solution (which there is, in this case it's simply what i'm generating the matrix from). When I get to use real data there will be no "true" solution.

As for which L2norm I want, I believe the difference lies in that matlab uses an L2 norm for a vector, whereas it uses a froebenius norm for a matrix. I believe that the latter is actaully the one I wanted. I will look into this PCG function.

Subject: fminsearch for functions

From: Matt J

Date: 21 Sep, 2011 21:06:28

Message: 10 of 12

"Karl " <karllandheer@gmail.com> wrote in message <j5dg7q$96v$1@newscl01ah.mathworks.com>...
>
> OK in my MLEM code I am just using the L2Norm as a kind of metric to show how close my reconstructed matrix is to my original matrix (I want it to be 0). I realize that my code will not, in general, find the minimal value of L2Norm but, correct me if i'm wrong, the most likely solution will be the one that minimizes the least squares IF there is a perfect solution (which there is, in this case it's simply what i'm generating the matrix from). When I get to use real data there will be no "true" solution.
================

Let me clarify. I don't know what MLEM even refers to when applied to a cost function like this. The term MLEM is usually reserved for Poisson loglikelihood minimization. Since I don't know what MLEM looks like here, I have no idea how it would behave, regardless of whether or not a perfect solution exists.

In any case, one other reason why "MLEM" might not be giving you the result you want, even assuming that a perfect solution exists and even assuming also that MLEM succeeds in minimizing the function, is that multiple minima may exist. If your process
A[p1,...,pN] is not a 1-1 transformation, or equivalently if the matrix Q is not full column rank, then norm(Q*radProf(:)-M(:)) will have a continuum of minima in addition to the simulated, perfect solution. Which solution your minimization algorithm converges to will depend on how you initialize.

Since we don't know how you circularly distribute your profile, that can't be assessed. Since your distributing circularly over a discrete rectangular grid, presumably it could involve some kind of rounding and/or nearest neighbour interpolation and this might introduce non 1-1 behavior. The only way to know for sure is to generate the Q matrix and try to compute its rank.

 
> As for which L2norm I want, I believe the difference lies in that matlab uses an L2 norm for a vector, whereas it uses a froebenius norm for a matrix. I believe that the latter is actaully the one I wanted. I will look into this PCG function.
=================

Then note that for a matrix, X

norm(X(:))=norm(X,'fro') %frobenius norm

but neither of these are equal to norm(X). So your code isn't doing what you intend.

Subject: fminsearch for functions

From: Matt J

Date: 21 Sep, 2011 21:13:28

Message: 11 of 12

"Karl " <karllandheer@gmail.com> wrote in message <j5dg7q$96v$1@newscl01ah.mathworks.com>...
>
> OK in my MLEM code I am just using the L2Norm as a kind of metric to show how close my reconstructed matrix is to my original matrix (I want it to be 0). I realize that my code will not, in general, find the minimal value of L2Norm but, correct me if i'm wrong, the most likely solution will be the one that minimizes the least squares IF there is a perfect solution (which there is, in this case it's simply what i'm generating the matrix from). When I get to use real data there will be no "true" solution.
================

Let me clarify. I don't know what MLEM even refers to when applied to a cost function like this. The term MLEM is usually reserved for Poisson loglikelihood minimization. Since I don't know what MLEM looks like here, I have no idea how it would behave, regardless of whether or not a perfect solution exists.

In any case, one other reason why "MLEM" might not be giving you the result you want, even assuming that a perfect solution exists and even assuming also that MLEM succeeds in minimizing the function, is that multiple minima may exist. If your process
A[p1,...,pN] is not a 1-1 transformation, or equivalently if the matrix Q is not full column rank, then norm(Q*radProf(:)-M(:)) will have a continuum of minima in addition to the simulated, perfect solution. Which solution your minimization algorithm converges to will depend on how you initialize.

Since we don't know how you circularly distribute your profile, that can't be assessed. Since your distributing circularly over a discrete rectangular grid, presumably it could involve some kind of rounding and/or nearest neighbour interpolation and this might introduce non 1-1 behavior. The only way to know for sure is to generate the Q matrix and try to compute its rank.

 
> As for which L2norm I want, I believe the difference lies in that matlab uses an L2 norm for a vector, whereas it uses a froebenius norm for a matrix. I believe that the latter is actaully the one I wanted. I will look into this PCG function.
=================

Then note that for a matrix, X

norm(X(:))=norm(X,'fro') %frobenius norm

but neither of these are equal to norm(X). So your code isn't doing what you intend.

Subject: fminsearch for functions

From: Karl

Date: 22 Sep, 2011 14:59:31

Message: 12 of 12

"Matt J" wrote in message <j5dk1o$npf$1@newscl01ah.mathworks.com>...
> "Karl " <karllandheer@gmail.com> wrote in message <j5dg7q$96v$1@newscl01ah.mathworks.com>...
> >
> > OK in my MLEM code I am just using the L2Norm as a kind of metric to show how close my reconstructed matrix is to my original matrix (I want it to be 0). I realize that my code will not, in general, find the minimal value of L2Norm but, correct me if i'm wrong, the most likely solution will be the one that minimizes the least squares IF there is a perfect solution (which there is, in this case it's simply what i'm generating the matrix from). When I get to use real data there will be no "true" solution.
> ================
>
> Let me clarify. I don't know what MLEM even refers to when applied to a cost function like this. The term MLEM is usually reserved for Poisson loglikelihood minimization. Since I don't know what MLEM looks like here, I have no idea how it would behave, regardless of whether or not a perfect solution exists.
>
> In any case, one other reason why "MLEM" might not be giving you the result you want, even assuming that a perfect solution exists and even assuming also that MLEM succeeds in minimizing the function, is that multiple minima may exist. If your process
> A[p1,...,pN] is not a 1-1 transformation, or equivalently if the matrix Q is not full column rank, then norm(Q*radProf(:)-M(:)) will have a continuum of minima in addition to the simulated, perfect solution. Which solution your minimization algorithm converges to will depend on how you initialize.
>
> Since we don't know how you circularly distribute your profile, that can't be assessed. Since your distributing circularly over a discrete rectangular grid, presumably it could involve some kind of rounding and/or nearest neighbour interpolation and this might introduce non 1-1 behavior. The only way to know for sure is to generate the Q matrix and try to compute its rank.
>
>
> > As for which L2norm I want, I believe the difference lies in that matlab uses an L2 norm for a vector, whereas it uses a froebenius norm for a matrix. I believe that the latter is actaully the one I wanted. I will look into this PCG function.
> =================
>
> Then note that for a matrix, X
>
> norm(X(:))=norm(X,'fro') %frobenius norm
>
> but neither of these are equal to norm(X). So your code isn't doing what you intend.

Let me say that I really appreciate your help and my linear algebra is not as good as it really should be (but this stuff is actually pretty interesting).

Anyways my question would be how to generate my matrix, then. I have a function which can determine Q*x, and I have Q*x = b, so in theory Q = b*inv(x), or, likely, Q = b*pinv(x), however keeping in mind that these are fairly large matrices I ran into a memory issue. Generally computing inverses is a no-no.

Additionally I am confused why calling pcg on my function does not work, it returns all 0s (which is not a solution). Again, thanks very much for your help.

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us