Solving a constrained system of linear equations

Question

I have a system of equations of the form y=Ax+b where y, x and b are n×1 vectors and A is a n×n (symmetric) matrix.

So here is the wrinkle. Not all of x is unknown. Certain rows of x are specified and the corresponding rows of y are unknown. Below is an example

| 10  |   |  5  -2  1 | | * |   | -1 |
|  *  | = | -2   2  0 | | 1 | + |  1 |
|  1  |   |  1   0  1 | | * |   |  2 |

where * designates unknown quantities.

I have built a solver for problems such as the above in Fortran, but I wanted to know if there is a decent robust solver out-there as part of Lapack or MLK for these types of problems?

---

My solver is based on a sorting matrix called pivot = [1,3,2] which rearranges the x and y vectors according to known and unknown

| 10  |   |  5  1 -2 | | * |   | -1 |
|  1  |   |  1  1  0 | | * | + |  2 |
|  *  |   | -2  0  2 | | 1 |   |  1 |

and the solving using a block matrix solution & LU decomposition

! solves a n×n system of equations where k values are known from the 'x' vector
function solve_linear_system(A,b,x_known,y_known,pivot,n,k) result(x)
use lu
integer(c_int),intent(in) :: n, k, pivot(n)
real(c_double),intent(in) :: A(n,n), b(n), x_known(k), y_known(n-k)
real(c_double) :: x(n), y(n), r(n-k), A1(n-k,n-k), A3(n-k,k), b1(n-k)
integer(c_int) :: i, j, u, code, d, indx(n-k)

    u = n-k
    !store known `x` and `y` values
    x(pivot(u+1:n)) = x_known
    y(pivot(1:u)) = y_known

    !define block matrices
    ! |y_known| = | A1  A3 | |    *    | + |b1|
    | |   *   | = | A3` A2 | | x_known |   |b2|

    A1 = A(pivot(1:u), pivot(1:u))
    A3 = A(pivot(1:u), pivot(u+1:n))
    b1 = b(pivot(1:u))

    !define new rhs vector
    r = y_known -matmul(A3, x_known)-b1

    % solve `A1*x=r` with LU decomposition from NR book for 'x'
    call ludcmp(A1,u,indx,d,code)
    call lubksb(A1,u,indx,r)

    % store unknown 'x' values (stored into 'r' by 'lubksb')
    x(pivot(1:u)) = r

end function

For the example above the solution is

    | 10.0 |        |  3.5 | 
y = | -4.0 |    x = |  1.0 |
    |  1.0 |        | -4.5 |

PS. The linear systems have typically n<=20 equations.

Answer 1

The problem with only unknowns is a linear least squares problem.

Your a-priori knowledge can be introduced with equality-constraints (fixing some variables), transforming it to an linear equality-constrained least squares problem.

There is indeed an algorithm within lapack solving the latter, called xGGLSE.

Here is some overview.

(It also seems, you need to multiply b with -1 in your case to be compatible with the definition)

Edit: On further inspection, i missed the unknowns within y. Ouch. This is bad.

Answer 2

First, i would rewrite your system into a AX=b form where A and b are known. In your example, and provided that i didn't make any mistakes, it would give :

    5 0 1     x1         13
A = 2 1 0 X = x2 and b =  3
    1 0 1     x3         -1

Then you can use plenty of methods coming from various libraries, like LAPACK or BLAS depending on the properties of your matrix A (positive-definite ,...). As a starting point, i would suggest a simple method with a direct inversion of the matrix A, especially if your matrix is small. There are also many iterative approach ( Jacobi, Gradients, Gauss seidel ...) that you can use for bigger cases.

Edit : An idea to solve it in 2 steps

First step : You can rewrite your system in 2 subsystem that have X and Y as unknows but dimension are equals to the numbers of unknowns in each vector. The first subsystem in X will be AX = b which can be solved by direct or iterative methods.

Second step : The second system in Y can be directly resolved once you know X cause Y will be expressed in the form Y = A'X + b'

I think this approach is more general.