Problem of NA in matrix

Hello, I need to do an exercice based on matrices on octave with the glpk function. Unfortunately octave returns me a matrix of NA surely due to a mistake in my code (there’s no NA in initial matrices). Here’s the code and some explanations:
The matrix M is a float matrix of size 100x1089
The matrix U is a float matrix of size 100x8
The matrix c represents the objective
The matrix rhs represents the right hand side of the constraints
ConstrCoef represents the left side of the constraints
The problem is surely in the lin_ls function
Here is the code:

1; % Script

% Vous devez implémenter ici votre méthode pour résoudre le problème
% 
%   min_{V >= 0}  M - U*V _1 tel que 

function V = votre_methodelinsn(M,U); 

[m,n] = size(M); 
[m,r] = size(U); 

% Résoudre le problème pour chaque pixel: pour tout j,
%   min{V(:,j) >= 0}  M(:,j) - U*V(:,j) _1
for j = 1 : n
    V(:,j) = lin_ls(M(:,j),U); 
endfor

endfunction


%  Résoudre les sous problèmes du type 
% 
%  min{x}  A*x - b _1 
% 
%  tel que  x >= 0. 

function x = lin_ls(b,A); 
 [m,n] = size(A); 
 c = [1; 1; 1; 1; 1; 1; 1; 1];
 s = ones(100,1);
 rhs = [b;s];
for j = 1 : n
  ConstrCoef(:,j) = [A(:,j);s]; 
endfor

% Formuler et résoudre le problème avec glpk
x = glpk(c,ConstrCoef,rhs);

endfunction

M = rand (100, 1089);
U = rand (100, 8);
V = votre_methodelinsn (M, U); 

Thank you for your question and example code.

I made a few modification to make your code runable on other systems as well and fixed one bug (calling the wrong function linls instead of lin_ls).

Can you provide a minimal and complete example (with matrices M and U, maybe as compressed mat-files). Something trivial like:

M = rand (100, 1089);
U = rand (100, 8);

Does not work to test the code :sweat:

Thanks a lot for your answer!
Here’s a file containing what you’ve asked for, hope it’ll help!
Hubble.mat (617.9 KB)

Thank you for the data.

Taking a closer look at the program, the intention is to solve 1089 linear programs (LPs), each of the form

min     x
s.t.  A*x  = b
        x >= 0

Is this right?

Can you explain why do you extend A by ones(100,100) (=ConstrCoef) and b by ones(100,1) ( =rhs)?

In general you a contructing an overdetermined LP, which most likely has no solution, even without the extension by ones. Reading the GLPK return value is always 10 (GLP_ENOPFS) No primal feasible solution:

[x,~,errnum] = glpk(c,ConstrCoef,rhs);

The returned NAs are the default return values in Octave’s GLPK wrapper. If GLPK computes a “useful” solution those NA get replaced by that solution.

This is just another way of telling the user: the desired computation failed.

To resolve the issue, it would be helpful to tell what exact problem you would like to solve with GLPK, especially the idea of the ones extension is not clear to me.

The form of the LP’s I want to solve is

min     x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8

    A(:,1)x1 + A(:,2)x2 + ... + A(:,8)*x8 = b
    x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 = 1

    x1, x2, x3, x4, x5, x6, x7, x8 >= 0

To call the glpk function, I need the objective function coefficients: c = [1; 1; 1; 1; 1; 1; 1; 1]

I also need the right-hand side value of each constaint: rhs = [b;1] (in the exercise, b is M(:,j) with j = 1 : 1089 )

And finally, I need the constraints coefficients.
So the matrix would be:

ConstrCoef =  A(:,1)  A(:,2)  A(:,3)  ... A(:,8)
                  1       1       1   ...     1 

(the A columns and the ones are in the same matrix)

To put the elements in the matrix, I put [A(:,j), 1] at the position ConstrCoef(:, j) with j = 1 : n (here, n = 8)

I don’t know how to update the code on the forum so I’ll repost it. I deleted s = ones(100,1) because I think it’s incorrect to use a matrix to only represent the value 1.
corrected code:

1; % Script

% Vous devez implémenter ici votre méthode pour résoudre le problème
% 
%   min_{V >= 0}  M - U*V _1 tel que 

function V = votre_methodelinsn(M,U); 

[m,n] = size(M); 
[m,r] = size(U); 

% Résoudre le problème pour chaque pixel: pour tout j,
%   min{V(:,j) >= 0}  M(:,j) - U*V(:,j) _1
for j = 1 : n
    V(:,j) = lin_ls(M(:,j),U); 
endfor

endfunction


%  Résoudre les sous problèmes du type 
% 
%  min{x}  A*x - b _1 
% 
%  tel que  x >= 0. 

function x = lin_ls(b,A); 
 [m,n] = size(A); 
 c = [1; 1; 1; 1; 1; 1; 1; 1];
 rhs = [b;1];
for j = 1 : n
  ConstrCoef(:,j) = [A(:,j);1]; 
endfor

% Formuler et résoudre le problème avec glpk
x = glpk(c,ConstrCoef,rhs);

endfunction

M = rand (100, 1089);
U = rand (100, 8);
V = votre_methodelinsn (M, U);

Thank you for the update :+1:

Below each post there is usually a small grey pencil icon “edit this post” as tooltip text. Or if it is hidden, see here How can I edit a post here on this discourse? - Site Feedback - Gitea.

The number of equality constraints is shrunk, but still the linear system matrix is very overdetermined [108 x 8] and my previous comment remains.

What makes you believe, that your system has a solution (paper, theorem, works in another software)? Can you relax the constraints to inequalities?

A*x <= b

Or does you model require equalities?

Hello!
First of all sorry for not responding for a few days.
I’ve had a chat with my teacher about this and in fact my approach of the problem is wrong (just as you said). So unfortunately i need to solve this on my own before using octave.
Thanks for helping the past few days!