3D and 2D matrix multiplication in the objective function

I have a 3D A matrix NxNxM and 2D B matrix (NxN), which will be an integer optimization variable matrix for the solver. I want to multiply A and B to obtain (NxNXM) matrix. Then I need to sum all of the elements of the resulting 3D matrix.

My objective function is the following:

ndesign = optimproblem('Objective',sum(A.*B,'all'));

I can do this multiplication if I don't use the solver. However, I receive the following error when I put it in the objective function:

Error using optim.internal.problemdef.ElementwiseOperator/checkIsValid
Argument dimensions 1454-by-1454-by-50 and 1454-by-1454 must agree.

Error in optim.internal.problemdef.Times/checkIsValid

Error in optim.internal.problemdef.Times.getTimesOperator

Error in .*

Error in networkdes (line 290)
ndesign = optimproblem('Objective',sum(A.*B,'all')); % Create optimization problem

I also need to convert my 3D matrix to sparse 3D matrix. I'm using ndSparse library to create it but I could not use it in the objective function again. Thanks in advance.

Edit for more explanations. Here is my optimization script:

A = %(1454x1454x50) 3D matrix 
A_sparse = ndSparse(A);
B = optimvar('b', [length(A),length(A)], 'Type', 'integer', 'LowerBound',0,'UpperBound',1);

ndesign = optimproblem('Objective',sum(A.*B,'all')); % Create optimization problem
% ndesign = optimproblem('Objective',sum(A_sparse.*B,'all')); % Create optimization problem for sparse A
ndesign.Constraints.const1 = ( B(diag(B)) == 0 );  % B_ii == 0
ndesign.Constraints.const2 = ( sum(B,'all') == 20 );

opts = optimoptions('intlinprog','Display','off','PlotFcn',@optimplotmilp);
[sol,fval,exitflag,output] = solve(ndesign,'options',opts);

[Stephan]

In R2019b and R2020a this works properly. Did you update your R2019a to the latest Update 9? Did you try to formulate the problem solver based?

[ebon T]

@Stephan Thanks for the advice Stephan. I updated my MATLAB to 2020b. However, I still get the same error.

[ebon T]

I edited the question with more details please have a look at it.

[Stephan]

The problem appears to occur only if there is the 3rd dimension in A:

A = rand(1454,1454,50);
B = optimvar('b', [length(A),length(A)], 'Type', 'integer', 'LowerBound',0,'UpperBound',1);
ndesign = optimproblem('Objective',sum(A.*B,'all'));

leads to error:

Error using optim.internal.problemdef.ElementwiseOperator/checkIsValid
Argument dimensions 1454-by-1454-by-50 and 1454-by-1454 must agree.

Error in optim.internal.problemdef.Times/checkIsValid

Error in optim.internal.problemdef.Times.getTimesOperator

Error in .*

Error in Untitled (line 4)
ndesign = optimproblem('Objective',sum(A.*B),'all');

Debugger does not work here, because we dont jump int the times function, when the error occurs. If you leave the 3rd dimension away, then it works without throwing an error.

A = rand(1454,1454);
B = optimvar('b', [length(A),length(A)], 'Type', 'integer', 'LowerBound',0,'UpperBound',1);
ndesign = optimproblem('Objective',sum(A.*B,'all'));

Also the normal calculations work as expected:

A =randi(10,10,10,3);
B =randi(10,10,10);
A.*B

gives:

ans(:,:,1) =

     9     2    30    45     8     4    12    18    18    25
     5    12     9    56     8     8     1    24    12    40
    81    35   100    21    80    10    18    27    18    30
    63    21     4    81    12    24     6     9     6    16
    30    20    40    42    54    16    63    36    25    80
    80    21     9    35    80    12    35     5    14    81
     9    63    54    30     9    50    42    80     7    15
     1    24    72    81    32    10    10    24    70    15
    42    42    24     2    10    48    90     7    30    20
     4     9    18     2    20    12    14    90    40    14


ans(:,:,2) =

    54     6    30    10    64    32     8    24    18    20
     1     3    30    24    16    64     3    80    18     5
    72    28    90    14    60     6    30    27    45    30
    21    56     3    45    24     9    12    27    30     8
    10   100    30    21    45    10    21    81    15    80
    80    63    15    63    48     4    49    45    16    72
    15    18    18    20    36    30    60    20    10    35
     7    54    56    36    28    15    20    60    42     9
    21    48    12     2     2    32    27    42    36    20
     3     5     6     4    45    30    20    90    36    21


ans(:,:,3) =

    72     2     6    35    64    32     4    21    18    50
     2     3     9    56     4    24     7    32    10    10
    18     7    30    35    30     4    21    81     9    10
    70    14     9    54    24     6    18    27    18    40
    45   100   100     7    81     2    56     9     5    80
    56    70    12    14     8    16    63    30    10     9
     3    45    45    20    36    10    54    10     2    15
     7    48    16    27     8    10    30    18    28    30
    49    36    30     8     8    40    27    49    30     4
     5    10    18     8    30    18    12    27    40    35

Since i have not much expirience with problem based approach, i dont know if this is a bug or if this is not possible. Maybe you can rewrite as solver based problem.

Maybe also Walter can help or even tell us if this is a bug, or whats going wrong here.

[ebon T]

Thank you for the effort Stephan. Yes it works with 2D matrix multiplication. I don't understand why it doesn't work with 3D matrix in the solver. I'll be gretaful if you can tag Walter to this question.

[Stephan]

I wrote a mail to Walter - im highly interested in his oppinion regarding this question. I think he will have a look at this soon. Maybe before he reads my mail...

[Walter Roberson]

No implicit expansion defined for that datatype.

[ebon T]

Is there a way to do this with different datatypes?

[Walter Roberson]

If you are using problem based optimization then that datatype is your only option. However if I understand correctly you can set up a for loop to construct the objective.

[ebon T]

Thanks a lot for the help !

I reshaped my A matrix to (n*m, n) and B is still (n,n )

A_re = reshape(A, [size(A,1)*size(A,3), size(A,1 )]); 

and changed my objective function to

ndesign = optimproblem('Objective',sum(A_re*B,'all ')); 

and it worked! However, I am not able to use ndSparse class sparse A matrix to do the same

A_sparse = ndSparse(A_re ); 
ndesign = optimproblem('Objective',sum(A_sparse*B,'all ')); 

When I run above script, I get the following error :

Error using double 
Conversion to double from optim.problemdef.OptimizationVariable is not possible.
Error in ndSparse>mkCompat (line 2803 ) 
X=double(X);
Error in * (line 1083 ) 
obj = finalObject( mkCompat(L)*mkCompat(R));
Error in networkdes_rev (line 20 ) 
ndesign = optimproblem('Objective',sum(A_sparse*B,'all ')); 

How can I do this and use it in the objective function?