How can I optimize the value of a matrix?

I want to find the optimal B matrix (2*4 matrix) that makes the elements of beta_d (1*4 vector) which is a function of B matrix, equal to the corresponding ones of a "given" beta_u (1*4 vector), i.e. I want to obtain the value of B matrix that makes beta_d(1,1) = beta_u(1,1) && beta_d(1,2) = beta_u(1,2) && beta_d(1,3) = beta_u(1,3) && beta_d(1,4) = beta_u(1,4).

On my attached code, as I am a beginner in MATLAB, I generated random value of B matrix and used 'while loop' to check if my condition is satisfied or not. This way is not efficient and takes very long time to give me a result.

Can anyone help me to optimize the value of B matrix in more efficient and faster way?

% Given H (4*2 matrix)
% Given A (2*2 diagonal matrix)
% Given C (4*4 matrix)
% Given P (4*4 diagonal matrix)
% Given C_l = 4; (constant value)
% Given beta_u (1*4 vector)

% I want to find the optimal B matrix (2*4 matrix) that makes the elements of the 1*4 beta_d vector equal to the corresponding ones of beta_u 1*4 vector
% i.e I want beta_d(1,1) = beta_u(1,1) && beta_d(1,2) = beta_u(1,2) && beta_d(1,3) = beta_u(1,3) && beta_d(1,4) = beta_u(1,4) 

%% My non-efficient code:

        B = zeros(2,4); % intial value of B 2*4 matrix
        beta_d = zeros(1,4); %intial value of beta_d 1*4 vector
        
        % compare each element of beta_d with the corresponding one in beta_u
        % optimize the value of B
        
        while(beta_d(1) ~= beta_u(1) || beta_d(2) ~= beta_u(2) || beta_d(3) ~= beta_u(3) || beta_d(4) ~= beta_u(4)) %check the condition
            
            B = randn(2,4); % generate a random 2*4 B matrix
            
            for j=1:1:2     % calculcate the 2*2 diagonal D matrix whose value depends on B 
                d(j) = (B(j,:)*P*B(j,:)')/((2^(2*C_l))-(norm(A(:,j))^2));
            end
            D = diag(d);   
            
            for i=1:1:4   % calculate beta_d (1*4 vector) whose value depends on B and D  
                V_d(i)=C(i,:)*P*B'*H(i,:)'*inv(1+H(i,:)*(A'*D*A+B*P*B')*H(i,:)');
                sigma_d(i)=norm((V_d(i)*H(i,:)*B-C(i,:))*(P^(1/2)))^2+(V_d(i)^2)*(1+H(i,:)*A'*D*A*H(i,:)');
                beta_d(i)=((P(i,i))/sigma_d(:,i)); 
            end
        end                

[Alexandra McClernon Ownbey]

you should initialize everything that is inside a for loop which means your variable will initialize for each iteration of the while loop. It also looks like you are checking your guess against a known value of integers. This also takes a long time since you don't give it a tolerance.

[Ahmed Nasr]

I set a tolerance and it takes a very long time too!

How can I use non-linear optimization to solve my problem?

Answer 1

[Rik]

The code below will fit a value for B with the fminsearch function. Note that it is relatively sensitive to a good initial guess. There are still several warnings mlint is giving you, but I will leave those for you to resolve.

edit: added the actual definitions for the input parameters

H = [-0.3208  -0.9784; -1.5994   -1.4689; -1.5197   -0.4568; -0.0993   -0.7667]; % 4*2 matrix
A = [-1 1; 0 1]; % 2*2 matrix
C = [-0.20  0.4  0.6 -0.2; -0.2  0.4  0.6 -0.2; 0.4  0.2 -0.2  0.4; 0.4  0.2 -0.2  0.4]; % 4*4 matrix
P = [250000   0   0   0; 0  250000  0   0; 0   0  250000  0; 0  0  0  250000]; % 4*4 matrix
C_l = 4; % constant value
beta_u = [50.2207   50.2207   20.3433   20.3433]; % 1*4 vector

%store inputs to a struct for shorter syntax
s=struct;[s.H,s.A,s.C,s.P,s.C_l]=deal(H,A,C,P,C_l);

initial_guess=randn(2,4);
OLS=@(b,input_vars)sum((myfun(b,input_vars) - beta_u).^2);%ordinary least squares cost function
opts = optimset('MaxFunEvals',50000, 'MaxIter',10000);
fit_val = fminsearch(OLS, initial_guess, opts,s);

function beta_d=myfun(B,input_vars)
%calculate beta_d from B and the other inputs

%load parameters
s=input_vars;[H,A,C,P,C_l]=deal(s.H,s.A,s.C,s.P,s.C_l);

beta_d = zeros(1,4); %intial value of beta_d 1*4 vector
for j=1:1:2     % calculcate the 2*2 diagonal D matrix whose value depends on B
    d(j) = (B(j,:)*P*B(j,:)')/((2^(2*C_l))-(norm(A(:,j))^2));
end
D = diag(d);
for i=1:1:4   % calculate beta_d (1*4 vector) whose value depends on B and D
    V_d(i)=C(i,:)*P*B'*H(i,:)'*inv(1+H(i,:)*(A'*D*A+B*P*B')*H(i,:)');
    sigma_d(i)=norm((V_d(i)*H(i,:)*B-C(i,:))*(P^(1/2)))^2+(V_d(i)^2)*(1+H(i,:)*A'*D*A*H(i,:)');
    beta_d(i)=((P(i,i))/sigma_d(:,i));
end
end

[Ahmed Nasr]

Thank you for your help.

As in the attached photo, this is not the result that I want. I want to obtain the value of B matrix that makes beta_d(1,1) = beta_u(1,1) && beta_d(1,2) = beta_u(1,2) && beta_d(1,3) = beta_u(1,3) && beta_d(1,4) = beta_u(1,4). If it's difficult to obtain equal values, it's acceptable for me to obtain slightly beta_d 's elements greater that beta_u 's ones. But, in the results, I obtained beta_d 's lower that beta_u 's ones. Can you help me in order to get the desired result?

[Rik]

If you don't share the input data as a mat file I can't really help you more. The only thing I can suggest is that you try other fitting functions.

[Ahmed Nasr]

This is my input data. If you need any further information, please let me know.

H = [-0.3208  -0.9784; -1.5994   -1.4689; -1.5197   -0.4568; -0.0993   -0.7667]; % 4*2 matrix
A = [-1 1; 0 1]; % 2*2 matrix
C = [-0.20  0.4  0.6 -0.2; -0.2  0.4  0.6 -0.2; 0.4  0.2 -0.2  0.4; 0.4  0.2 -0.2  0.4] % 4*4 matrix
P = [250000   0   0   0; 0  250000  0   0; 0   0  250000  0; 0  0  0  250000] % 4*4 matrix
C_l = 4; % constant value
beta_u = [50.2207   50.2207   20.3433   20.3433]; % 1*4 vector

I want to obtain the value of B matrix that makes beta_d(1,1) = beta_u(1,1) && beta_d(1,2) = beta_u(1,2) && beta_d(1,3) = beta_u(1,3) && beta_d(1,4) = beta_u(1,4). If it's difficult to obtain equal values, it's acceptable for me to obtain slightly beta_d 's elements greater that beta_u 's ones. But, in the results using the above answer, I always obtain beta_d 's lower that beta_u 's ones. Can you help me in order to get the desired result?

[Rik]

Do you have a B that returns a better fit? Because your original code generates a new random matrix every iteration it is taking ages to find a suitable one, even if I use while sum(abs(beta_d-beta_u))>=1

If you have a good B I can more easily experiment with what fitting functions would work.

[Ahmed Nasr]

No, I don't have. All what I know, B is a 2*4 matrix and I want to find its optimal value that makes "beta_d" vector equal to "beta_u" vector.

[Rik]

Are you certain it is possible? With the function you wrote it might be impossible to find a B such that the beta vector will be what you need it to be. I'm not sure there is an easy way to check. I do know fminsearch is relatively sensitive to your initial guess, as it can somewhat easily get stuck in a local minimum. If you have the curve fitting toolbox you will have more options. Note that any fitting algorithm may have issues fitting 8 free parameters.

What you could do is trying to adapt the cost function so it doesn't try to minimize the sum of squared difference, but looks at the relative offset. I doubt that would make a large difference, but it is something you could try.

[Ahmed Nasr]

We can use a constraint that:

Will this constrait help us to find the optimal B?

[Rik]

Is there really nothing in your background information that gives a hint as to what would be a valid intial guess? Was there a specific reason you were using randn at first? And are you sure there even exists a B that would result in a better fit? Do you have actual reasons for such a belief?

You can use the setup below to add the constraint to your fitting function. It will add inf to your cost function if the test_trace_condition function returns false. This results in more or less the same quality of fit.

H = [-0.3208  -0.9784; -1.5994   -1.4689; -1.5197   -0.4568; -0.0993   -0.7667]; % 4*2 matrix
A = [-1 1; 0 1]; % 2*2 matrix
C = [-0.20  0.4  0.6 -0.2; -0.2  0.4  0.6 -0.2; 0.4  0.2 -0.2  0.4; 0.4  0.2 -0.2  0.4]; % 4*4 matrix
P = [250000   0   0   0; 0  250000  0   0; 0   0  250000  0; 0  0  0  250000]; % 4*4 matrix
C_l = 4; % constant value
beta_u = [50.2207   50.2207   20.3433   20.3433]; % 1*4 vector

%store inputs to a struct for shorter syntax
s=struct;[s.H,s.A,s.C,s.P,s.C_l]=deal(H,A,C,P,C_l);

initial_guess=[-1 -1 1 1;1 -1 1 1]/2;%trial and error
if ~test_trace_condition(initial_guess,s)
    error('change initial condition to a valid value')
end

OLS=@(b,input_vars)sum((myfun(b,input_vars) - beta_u).^2);%ordinary least squares cost function
t_inf_f_0=@(tf)inf^(tf)-1;%convert true to inf and false to 0
costfun=@(b,input_vars) OLS(b,input_vars) ...
    + t_inf_f_0(~test_trace_condition(b,input_vars));
opts = optimset('MaxFunEvals',50000, 'MaxIter',10000);
fit_val = fminsearch(costfun, initial_guess, opts,s);

if test_trace_condition(fit_val,s) && ~isinf(costfun(fit_val,s))
    clc
    disp('fit succesful')
    beta_u %#ok<NOPTS>
    beta_d=myfun(fit_val,s) %#ok<NOPTS>
else
    clc
    disp('fit failed')
end

function beta_d=myfun(B,input_vars)
%unchanged
end
function tf=test_trace_condition(B,input_vars)
%returns true for a valid B
s=input_vars;[A,P,C_l]=deal(s.A,s.P,s.C_l);
d=zeros(2,1);
for j=1:1:2     % calculcate the 2*2 diagonal D matrix whose value depends on B
    d(j) = (B(j,:)*P*B(j,:)')/((2^(2*C_l))-(norm(A(:,j))^2));
end
D = diag(d);

tf= trace(B*P*B' + A*D*A') <= trace(P);
end

[Ahmed Nasr]

At first, I am really grateful for your help. According to your questions, Yes, I am sure that beta_d must equal beta_u at the end. One note that can help us, in order to obtain beta_u, I calculate B_u (4*2 matrix) using MMSE. Can we use intial guess = B_u'; ?

I tried your code and this is the result:

The problem in the last element of beta_d(1,4).

[Rik]

You keep forgetting: you have more data than I do. If you don't share the values of the variables you're using I can't help you at all. Either post the code that will generate B_u, or post the values.

As for your second question: did you try to determine which values in B have an influence on the fourth value of beta_d? I wouldn't be surprised if all have an influence. Maybe one of the values in B has a larger inpact than the others. You could try manually adjusting the values of B to see what effect they have. You could try reasonable ranges for all 8 values to see what would bring your output closer to the output.

[Ahmed Nasr]

I tried to generate random H (4*2 matrix). sometimes, the code works well and give me successful fit and on other random values of H, it gives a failed fit. do you know a reason for that?

[Ahmed Nasr]

I want to save only the results of the successful fit, How can I do that?

[Rik]

You know the values of your variables. You can save your variables as normal.

So if I understand you correctly, my answer solved your question now?

[Ahmed Nasr]

Yes, thank you for your efforts. I really apperiate it.

[Ahmed Nasr]

If I can ask you another question:

I have 4 users, and I want to simulate the outage probability which is the probability of the user data rate is less than the predefined threshold value i.e. if the threshold rate per user is = 0.512.

After changing the signal to noise ratio SNR from -20 dB to 30 dB, I obtain the rate matrix with dimensions (number of users × the length of SNR vector), where each element in this matrix represents the rate of certain user at specific SNR value as follows:

SNR_vector_dB=-20:10:30; %the SNR vector
R_th = 0.512; % the user data rate threshold value 
% Rate matrix with dimensions (number of users =4 * the length of SNR vector=6)
R = [0.0021    0.0191    0.1466    0.7368    1.5617    1.9838; ...
0.0021    0.0191    0.1466    0.7368    1.5617    1.9838; ...
0.0041    0.0370    0.3190    0.9267    1.6591    2.1856; ...
0.0041    0.0370    0.3190    0.9267    1.6591    2.1856];

Can you help me in simulating the outage probability versus SNR?