trying to find a linear combination of matrixes in order to minimize the error between the linear combination and a target matrix

Hi everyone, i'm trying to find a solution to my problem, i need to calculate coefficient for a linear combination of matrix in order to minimize the difference between a target matrix and my linear combination.

in symbols : Error = TargetMatrix - (aA+bB+cC+dD+eE+fF+gG+hH+iI+jJ+kK+lL);

I need to find a, b, c , d ,e , f ,g, h , i ,j ,k and l in order to obtain the minimum Error;

[Jon]

To minimize something you must say how you are measuring how big it is. In your case you are trying to minimize, a matrix, so you must say, what norm you are using to measure how "big" the matrix is. See https://www.mathworks.com/help/matlab/ref/norm.html

Answer 1

[John D'Errico]

Trivially easy. What you need to do it to convert the problem to a set of VECTORS. Unroll the matrices into vectors. Then make one large matrix. Next, you need to learn to use matrices in general. DON'T have 9 different matrices. Have ONE three dimensional matrix. LEARN TO USE MATLAB AS IT IS DESiGNED TO BE USED.

matrices = rand(4,4,9);

Lacking any data, I have created 9 random matrices, each 4x4.

M = reshape(matrices,[4*4,9]);

I'll choose an array of all ones as the goal matrix, so we can see how well we did.

Target = ones(4,4);
Treshape = Target(:);

Now I have ONE matrix, of size 16x9. Here, each of the original matrices is one column of M. The right hand side for the problem is just the target matrix, also reshaped into a vector of length 16,

At this point, the problem is just one of a simple multiple linear regression.

coeffs = M\Treshape

coeffs = 9×1

0.4437 0.3894 0.1487 0.2571 0.2808 -0.1499 -0.1652 0.2343 0.4149

How well did we do?

pred = reshape(M*coeffs,[4,4])

pred = 4×4

1.0693 1.0065 0.7762 1.0667 0.8101 1.0037 0.9704 0.8841 1.2350 0.8889 1.1612 0.9441 1.0293 1.1588 0.6914 0.9760

Not too bad, actually, since the target was all ones.

Again, learn to use MATLAB, as it was designed to be used.

[Peter Marfin]

Thank you so much, i had the idea of doing something like this but i'm quite new to matlab so i didn't realize it was that simple

[Jon]

@John D'Errico nice use of reshape and \, What matrix norm of the error=pred - Target does this approach minimize? I think it is the Frobenius norm, is that correct? The reason I am asking is I was coding up an alternative approach, that was more general (in terms of the norm used), but definitely a lot more elaborate.

[John D'Errico]

Yes, it would be the Frobenious norm that is minimized. Frobenious just computes the sum of squares of elements of the difference, and this approach would minimize that. There was no statement of how the error should be defined in the question, so I could only make an assumption. Anyway, I've found it to be pretty rare that people actually state a preference, and on a somewhat vague question where the only goal is to minimize error, I have a funny feeling that the implied norm minimized is the last thing to be in question.

If I wanted to solve the problem for other norms however, I would NOT just throw an optimizer at the problem. Instead, I would recall there are classic solutions for the 1-norm and infinity norm problems, just using linprog and slack variables. So no need to use a general optimizer, which then requires things like starting values. (Ok, I can probably write the 1-norm case from memory, and I'm pretty sure the inf-norm has a solution via linprog too, but it has been 40 years since I learned the trick.)

[Jon]

Wow, thank you so much for your insights on that! Good to know that for the other norms there are linear programming solutions. I was already starting to see that solving this with fminunc was getting a little fiddly for larger collections of matrices. Woud be much better to have a dependable linear programming solution. Regarding the choice of norm, I'm wondering if there is some notion that if we make one norm "very" small we also make the others small, so for most practical purposes it doesn't really matter what norm we choose. Maybe the choice of norm really only comes in to play for situations where some constraints keep the minimal norm relatively large, so the differences between norms become more apparent. Anyhow, I guess this discussion is probably going well beyond what was asked by the OP, but it is interesting, at least to me :)

[Bruno Luong]

@John D'Errico " I can probably write the 1-norm case from memory, and I'm pretty sure the inf-norm has a solution via linprog too, but it has been 40 years since I learned the trick.)"

For matrix actually the Inf-norm is the 1-norm of its transpose.

A=rand(4)

A = 4×4

0.1831 0.5361 0.3770 0.1296 0.9693 0.5670 0.1173 0.3901 0.3642 0.4925 0.6051 0.3450 0.6674 0.0536 0.2527 0.2504

norm(A,Inf)

ans = 2.0438

norm(A.',1)

ans = 2.0438

You probably have in mind vector where 1-norm and Inf-norm are fundamentally different.

Never minimize 1-norm or Inf-norm for matrix myself, but yes it should be doable with linprog as they are "flat" norm, and I wonder if problem-based solver can be smart enough to figure it out.

Answer 2

[Jon]

Here's a much more elaborate approach that can be modified to use a specific choice of matrix norm to be minimized. I would use @John D'Errico's very clean and simple approach unless there were some overwhelming reason that the choice of matrix norm was very important. I am including this just in case it is of interest, and also to get @John D'Errico perspective on whether there are ever any cases where some other matrix norm might be needed. I also not that this approach requires the Optimization toolbox. Also using fminunc here seems to have some potential issues when the number of matrices to be regressed gets large (e.g. 10).

% Find linear combination of matrices to approximate a target matrix
% min(||A - a1*B1+a2*B2+ ...anBn||) where a1,a2,.. are scalars, A, B1, B2,...
% are mRow by nCol matrices and || || is a matrix norm
% in this case we will use the p=2 norm for || ||, but you could edit this
% to use the norm of your choosing

% Make some example data using random matrices, you would use the matrices
% of interest for your problem
mRow = 3;
nCol = 3;
n = 7; % number of matrices to be included in fit

% create an example where we know the exact answer
B = randn(mRow,nCol,n);
bExact = randn(n,1);
Aexact = lincomb(B,bExact);
% Create target matrix that can't be exactly acheived as linear combination
A = Aexact + 100*eps*randn(mRow,nCol);

% Find optimal solution, you should check exit flag, for this example I
% skipped this for the example
[bOpt,exitflag] = bestLinComb(A,B)

                                                        First-order 
 Iteration  Func-count       f(x)        Step-size       optimality
     0           8          8.66188                          2.78
     1          16          2.69719        0.36031           3.08  
     2          32          1.46508       0.338786           2.27  
     3          48          1.21195       0.265759           2.63  
     4          64         0.721338       0.482196           1.32  
     5          80          0.58381       0.245134           1.24  
     6         104         0.548785       0.080194           1.96  
     7         120         0.226155       0.663275           1.18  
     8         136         0.198837       0.278846           2.41  
     9         152         0.103643       0.244825           1.63  
    10         168          0.08529       0.167117           1.59  
    11         184        0.0481025       0.353594           1.76  
    12         208        0.0292037      0.0550715            1.2  
    13         224        0.0196182       0.191943           1.95  
    14         248        0.0152428      0.0308406           1.09  
    15         264          0.01388       0.139991            2.7  
    16         288        0.0137119         0.0172           2.36  
    17         304         0.010171       0.176021           1.25  
    18         320       0.00511634       0.266131           1.31  
    19         336       0.00477305       0.151178           2.41  
                                                        First-order 
 Iteration  Func-count       f(x)        Step-size       optimality
    20         352       0.00329684       0.383929           1.93  
    21         376       0.00247798      0.0657629           1.35  
    22         392       0.00140657       0.200203           3.01  
    23         408       0.00124177       0.146318           1.45  
    24         424       0.00116343       0.117202           1.47  
    25         440      0.000939582       0.165542           1.83  
    26         464      0.000822947      0.0507037           2.09  
    27         480      0.000469699       0.239479           3.02  
    28         496      0.000269173       0.201411           1.77  
    29         520      0.000244038      0.0258355           1.44  
    30         536      0.000200083       0.214226           1.58  
    31         552      0.000120895       0.221581           2.96  
    32         568      8.95311e-05       0.212502           1.94  
    33         584      6.93162e-05       0.173591           1.96  
    34         600      6.30294e-05       0.143538            2.3  
    35         616      4.48097e-05       0.182732           2.94  
    36         632      2.83974e-05       0.419394           1.75  
    37         656      2.05193e-05      0.0265209           1.89  
    38         672       1.5644e-05       0.171718           2.83  
    39         688      1.19023e-05       0.234888           2.14  
                                                        First-order 
 Iteration  Func-count       f(x)        Step-size       optimality
    40         704      7.35839e-06       0.158638           1.46  
    41         728      6.62833e-06      0.0292406           1.34  

Optimization stopped because the norm of the current step, 4.611555e-07, is
less than options.StepTolerance = 1.000000e-06.

bOpt = 7×1

1.0403 0.2393 0.0296 -0.0579 0.3491 0.9840 -1.5403

exitflag = 2

% Compare optimal b and bExact
deltab = bOpt - bExact

deltab = 7×1

1.0e-05 * -0.0586 -0.1907 -0.2597 -0.0862 -0.2076 0.0274 -0.0580

% Compare approximated A and target
Aopt = lincomb(B,bOpt);
delta = A - Aopt

delta = 3×3

1.0e-05 * -0.1317 0.1742 0.2771 -0.0563 0.1065 0.4667 0.5254 -0.0488 -0.0866

% Compare to very simple approach
bAlt = reshape(B,mRow*nCol,n)\A(:)

bAlt = 7×1

1.0403 0.2393 0.0296 -0.0579 0.3491 0.9840 -1.5403

%
function [b,exitflag] = bestLinComb(A,B)
% A - target matrix
% B - m by n by p array of matrices to be combined
b0 = zeros(size(B,3),1);
options = optimoptions('fminunc','Display','iter-detailed',...
    'MaxFunctionEvaluations',1e6);
[b,~,exitflag] = fminunc(@errnorm,b0,options);
% function to find best linear combination of matrices
    function e = errnorm(b)
        % error whose norm is to be minimized
        % you could use a different norm here, e,g, 1, inf,"fro"
        e = norm(A - lincomb(B,b),2); % minimize svd
    end
end

function B = lincomb(A,a)
% calculates linear combination of matrices
% A is a m by n by p array of matrices
% a is a length p vector
% B = a(1)*A(:,:,1) + a(2)*A(:,:,2) + ... a(p)*A(:,:,p)

% compute linear combination, we first turn each matrix into a  vector
% columnwise, form a matrix from these vectors and then use ordinary matrix
% multiplication, and finally reshape to the original size
[m,n,p] = size(A);
B = reshape(reshape(A,m*n,p)*a(:),m,n);
end