How to write a better "if condition" to eliminate some unwanted matrices?

To be clear, I'm changing my question as follows with more details:

------

Here is the code segment I'm working on:

NphaseSteps = 4; 
% Assume that TM is a 10x10 random matrix that includes both real and complex numbers. 
% But it is fixed after created randomly.
TMsize = length(TM);
focusX = 5; focusY = 5;

function inFieldA = WSfunction(NphaseSteps, TM, TMsize, focusX, focusY)

% This function is supposed to return inFieldA that gives 1) the maximum value
% at the target position (focusX,focusY) inside the outField matrix 2) the 
% optimized outField such that all other elements other than (focusX,focusY) 
% inside outField have "at most" half of the target position value.

tic

fprintf('The number of possible random matrices: %s \n',NphaseSteps^(TMsize^2));

phases = exp( 2*pi*1i * (0:(NphaseSteps-1))/NphaseSteps );
nIterations = 5;
inField = cell(1,nIterations);

i = 1; j = 0;

while i <= nIterations %number of iterations

     ind = randi([1 NphaseSteps],TMsize,TMsize);
     inField{i} = phases(ind);

     j = j+1;

     outField{i} = TM * inField{i};
     outI = abs(outField{i}).^2;
     targetIafter(i) = abs(outField{i}(focusX,focusY)).^2;
     middleI = targetIafter(i) / 2;

     if max(max(outI)) == targetIafter(i)
        outI(focusX,focusY)=0;
           if (outI > middleI)  == zeros(TMsize,TMsize)
              i = i + 1;
           end
     end

     if mod(j,10^6) == 0
        fprintf('The number of random matrices tried: %d million \n',j/10^6)
     end         
end

fprintf('The number of random matrices tried in total: %d \n',j)

[maxInt, maxInd] = max(targetIafter);
inFieldA = inField(maxInd);

toc

end

This program basically creates a 10×10 random input matrix inField that only includes 4 elements (number of phase steps) and multiplies it with a 10×10 fixed matrix TM to obtain a 10×10 output matrix outField. Now, I want to fix two things about the program:

1. Since the inField matrices are created randomly, the program will create equal matrices many times. I want to eliminate the equal inField matrices before they go into the "if statement". For this I need to change the cellfun part but I don't know how to do it in the most efficient way.
2. In the outField matrix I want a specific matrix element to have the largest value. First part of the "if statement" deals with that. But, I additionally want the other elements (other than my target matrix element) to have values smaller than the half of the target value. In the second part of the "if statement" I wanted to handle this situation but I think it's too complicated and not very efficient. Can you offer a better way to do it?

[dpb]

"I want its maximum element to be at a desired target location say (5,5)."

What does this mean, precisely?

Set outI(targetIafter) to whatever value max(outI(:)) is?

Or keep outI(targetIafter) as found but halve everything else?

Or something else?

A small illustrative sample would be informative.

[Sahin Kurekci]

@dpb I've edited my question by giving more details. I hope my intention is clear now.

[dpb]

Not really enough to code, no (at least for me, maybe someone else can follow, I dunno 'bout that).

As suggested before an illustrative same would likely help.

[Bruno Luong]

Since the inField matrices are created randomly, the program will create equal matrices many times. I want to eliminate the equal inField matrices before they go into the "if statement". For this I need to change the cellfun part but I don't know how to do it in the most efficient way.

Why make you think of that? The number of possible combinations is

6^(10*10) > 6.5e77

you think with 5 iterations you will have a chance to get 2 exact same matrix? I think winning lotery is easier.

You lost big percent of CPU time to take of case that almost never happens.

[Sahin Kurekci]

@dpb I've edited my code again for more clarity. Now, you only need a random complex matrix TM to work on it. You can see the results after graphing the outField matrix for the optimized inField that is given by the program. @Bruno Luong, you are right! But, initially my program was taking too much time to execute and I thought that since I have two conditions inside "if statement", equal inField matrices are causing that trouble while they are inside "if".

[Bruno Luong]

No what take time is that your code try of growing the cell, making a storage and then comparing, all that effort would never detect repeat generate matrix (the probability is tiny, tiny, tiny...). That's my point.

Focus on the rest rather than bother with unique creation of random matrix.

[Bruno Luong]

Try to remove both cell storages, and see whereas your code run faster. Still if you don't see solution is cached for many iterations, your criteria is probably too strict (that is what I suspect the case).

[Sahin Kurekci]

@Bruno Luong, I also suspected of my criteria being strict but to get rid of this I changed the variable NphaseSteps to many other values. In theory increasing its value should cause variety inside the inField matrix and so it will be easier to catch a maximum value at the desired target point in outField. (This is my opinion at least.) But, it also took so much time and at one point I ended the execution.

I actually tried to give a short, clear definition of my problem in another question, you can check the following link and I hope the conditions I want to implement inside the program will be clearer to you: How to find an optimized matrix by changing the elements of another matrix?

[Bruno Luong]

I don't know about your problem, but multiply by a matrix is usually a regulator operator (integral), unless the matrix is very special. Usually the product will not be a sharp peak at one target point.

Here is an idea: You can always solve for continuous matrix (INFIELD)

TM * INFIELD = M_target

where M_target is defined as M_target(5,5) = 1, and 0 elsewhere.

Meaning for unconstraint case

INFIELD = TM \ M_target

Or introduce that fact that abs(phase)=1; using fmincon to minimize

f(INFIELD ) := norm(TM * INFIELD - M_target,'fro')^2
under constraint 
abs(INFIELD_(i,j)) are all equal, for all i,j from 1 to 10.

Then try to see which of you discrete phase matrix solution meets your criteria.

[Sahin Kurekci]

@Bruno Luong, thank you for your suggestion. I actually tried your method of dividing TM to outField but abandoned this idea immediately since there was a very low possibility of finding an inField matrix which is formed from the exact 4 numbers I want. But, looking at the correlations is a very nice idea. I can try that, thank you again.

[Bruno Luong]

I just edit my comment with the idea of minimization with constraint.

[Sahin Kurekci]

@Bruno Luong can you expand your answer further? I constructed the outField matrix with target location value being equal to 1 and all the others are 0, then I did the division. Now, should I write a new function that finds the norm of the difference between TM*inField and outField? Is the last line you wrote is supposed to be my constraint for this function? And lastly, where will I use fmincon here? Sorry for all these questions but I'm in a hurry to finish this program and I'm also new to MATLAB programming; all these make things complicated for me.

Answer 1

[Bruno Luong]

Sahin: there is one thing obvious I didn't see until now, but the each column of the product

TM * inField

depends independenly of column of inField (since it matrix product), there for if you want to have a goal of outField(5,5) to stick out, all you need to do is solving two subproblems of (matrix x vectors) form.

TM*inGoal is proportional to e5:=[0,...0,1,0...]' (1 at 5th position)
TM*outgoal is as small as possible

where inGoal, outGoal are both (10x1) vectors.

An approximation of inGoal is very to find, all you do is small subsub-problem:

TM(:,j)*ingoal(j) maximizes a correlaton with e5

For outGoal, what you need is to minimize

argmin|TM*x|_inf / |x|_inf, where |x|_inf := max(abs(x(i)).

If we denote y:=TM*x, this is equivalent to minimizing

argmin |y|_inf / |inv[TM]*y|_inf

which is the same as maximizing

argmax |inv[TM]*y|_inf / |y|_inf.

This is nothing than find the matrix norm inv[TM]_inf

The solution of which is (just look at google):

argmax |inv[TM]*y|inf / |y|_inf = max_i sum_j (inv[TM]_i,j)

and the proof of that also provides the solution y that make the ratio reach the matrix norm.

One you have y, just compute x = TM \ y, and that gives pretty much the phase vector outGoal to make TM*outGoal small.

Put together to get

inField = [outGoal, ..., outGoal, inGoal, outGoal, ..., outGoal]

where the inGoal is at 5th column, there you are.

[Bruno Luong]

Following Matt's idea where I modify to 4th power, and proper normalize, and introduce the parametrization

phase = atan2(2t/(1+t^2),(1-t^2)/(1+t^2))

here is the code using fminunc (and not need for fmincon) to find the phase (exp(i*phase) is the inGoal vector) the makes one component of (TM*in) dominant:

 n = 10;
 rowgoal = 5;
 TM = randn(n,n)+1i*rand(n,n);

% Try to estimate p such that, zin = exp(1i*p) such that inres := abs(Tm*zin)
% has the element inres(rowgoal) dominant

t0 = zeros(n,1);
t = fminunc(@(t) funin(t,TM,rowgoal), t0);
zin = t2phase(t);
inres = abs(TM * zin)

% Try to estimate p such that, zout = exp(1i*p) such that outres := abs(Tm*zout)
% has small elements

t0 = zeros(n,1);
t = fminunc(@(t) funout(t,TM), t0);
zout = t2phase(t);
outres = abs(TM * zout)

cla
plot(inres,'b')
hold on
plot(outres,'g')
plot(rowgoal+[0 0],ylim(gca),'-.r');

function d = funin(t,TM,rowgoal)
n = size(t,1);
z = t2phase(t);
y = TM * z;
y4 = (y .* conj(y)).^2;
d = (n-1)*(y4(rowgoal));
y4(rowgoal) = 0;
d = sum(y4)-d;
end

function d = funout(t,TM)
n = size(t,1);
z = t2phase(t);
y = TM * z;
y4 = (y .* conj(y)).^2;
d = sum(y4);
end

function z = t2phase(t)
t2 = t.^2;
s = 1+t2;
sn = 2*t./s;
cs = (1-t2)./s;
z = cs + 1i*sn;
end

You can observe that for random TM matrix, it's very difficult to make the 5th component of the output twice larger than the others.

I don't know if there is local minima that prevents minimization algorithm to converge to the best solution, you might try to run fminunc in multi-starting mode and peak the best solution.