How to solve a linear least squares problem in a dense/sparse format

New to AMPL
1

Hello! I want to know the general AMPL code to solve a linear least squares problem.
I think it’s quite tricky for me to write the AMPL code for this problem.

Here is the assumption and the equation:

A \in \mathbb{R}^{m \times n } \\ \mathbf{x} \in \mathbb{R}^{n \times k} \\ \mathbf{b} \in \mathbb{R}^{m \times k} \\ \underset{\mathbf{x}}{\min} || A\mathbf{x} - \mathbf{b} ||^2

If we have this data declaration, I think the matrices are in the dense format.

set M;
set N;
set K;
param A {M, N};
param b {M, K};
var X {N, K};

The question is how I can define my minimization function for the problem? Is there any way to use functions that will work for matrix or vector representation? If not, then can I store the intermediate results for the minimization?

My question is actually to ask how I can get each row of A\mathbf{x} to subtract it by \mathbf{b}. I want to do like this:

minimize length :       
sum{i in M} # sum of each vector for each row
   
   {j in N, k in K} # just for indexing 
   row_Ax := A[i, j] * X[j, k] # store the row i of Ax for the matrix multiplication.

  sum {k in K} # to get the length of (Ax - b)
  (row_Ax[k] -  b[i, k]) * (row_Ax[k] -  b[i, k])  # (the actual minimization code)

I think it’s not the way of the AMPL. But can anyone let me know how to solve this and even solve this one in a sparse matrix format? I don’t know it’s possible to declare a sparse matrix and do the operations to solve the linear least squares problem with the AMPL.

1 Like
2

Hi @lch32111,

It should not be really tricky to write this problem into ampl after you do the computations.

To define the minimization function let’s first compute the elements for the matrix Ax-b, which is a mxk-dimensional matrix. The values in Ax are:

(Ax)_{i,j} = \sum \limits_{c=1}^{N} a_{i,c} \cdot x_{c,j}

Then, values for Ax-b are:

(Ax-b)_{i,j} = \sum \limits_{c=1}^{N} a_{i,c} \cdot x_{c,j} - b_{i,j}

The translation into ampl is straightforward (following your model notation):

    sum{c in N} A[i,c]*X[c,j] - b[i,j]

Now you can use these elements to write the minimization function. I am not sure which matrix norm are you using, for example, if you want to minimize the sum of the squared differences over the elements (squared euclidean norm), it would be something like:

minimize sum_sq_element:
    sum{i in M, j in K} (sum{c in N} A[i,c]*X[c,j] - b[i,j])^2;

Ampl stores your matrix in a sparse format, so you do not have to worry about that. You may specify the matrices A and b in a data file, or maybe use some dictionary in some programming language + ampl api.

In the following model, non specified values for A and b are assumed to be zero.

param m;
param n;
param k;
set M := 1 .. m;
set N := 1 .. n;
set K := 1 .. k;

param A {M, N} default 0;
param b {M, K} default 0;
var X {N, K} >= 0;

minimize sum_sq_element:
    sum{i in M, j in K} (sum{c in N} A[i,c]*X[c,j] - b[i,j])^2;

Possible data file.

data;

param m := 3;
param n := 2;
param k := 3;

param A :=
1 1 1
1 2 3
2 2 1.5
3 1 -1
;

param b :=
1 1 1.5
1 2 2
1 3 1
2 1 -1.5
2 2 1
2 3 -1
3 1 0
3 2 10
3 3 3.5
;

If you are using an api, you can just specify the data like (Python example):

A={(1,1):1, (1,2):3, (2,2):1.5, (3,1):-1}
ampl.param['A'] = A

Similar for the b matrix.

Then pick your solver and display the solution:

option solver gurobi;
solve;
display X;

Hope it helps!

1 Like
3

I didn’t know that way to solve it.

You answer is so perfect enough to make me know what I will do.

I appreciate your help @marcos !

1 Like