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.

the usual approach is to perform the required algebraic manipulation to put the system into standard form, then use a standard solver. — agentp
@agentp I think this is what I show above. I re-arrange the equations in block form in order to bring it into Ax=b form. But, is there common way to do this. To me this is general enough problem to warrant more formalism. — John Alexiou

sascha sascha · Accepted Answer · 2017-04-17T12:46:07

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.

Solving a constrained system of linear equations

2 Answers