求解受约束的线性方程组

Question

我有一个形式的方程组y=Ax+b，其中y，x和b是n×1向量，A是n×n(对称)矩阵。

所以这是皱纹。并非所有x都是未知的。指定了x的某些行，而y的相应行是未知的。下面是一个示例

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

其中*指定未知数量。

我已经在Fortran中为上述问题构建了一个求解器，但我想知道是否有一个像样的健壮求解器作为Lapack或MLK的一部分来解决这些类型的问题？

---

我的求解器基于称为pivot = [1,3,2]的排序矩阵，该矩阵根据已知和未知重新排列x和y向量

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

以及使用块矩阵解决方案和LU分解求解

! 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

对于上面的示例，解决方案是

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

附言。线性系统通常具有n<=20方程。

Answer 1

只有未知数的问题是一个线性最小二乘问题。

您的先验知识可以通过相等约束(修复某些变量)引入，将其转换为线性相等约束的最小二乘问题。

lapack中确实有一种算法可以解决后者，称为xGGLSE的算法。

以下是一些概述。

(看起来，在你的情况下，你需要将 b 乘以 -1 才能与定义兼容)

编辑：在进一步检查时，我错过了y内的未知数。哎哟。这很糟糕。

Answer 2

首先，我会将您的系统重写为 AX=b 形式，其中 A 和 b 是已知的。在您的示例中，如果我没有犯任何错误，它将给出：

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

然后，您可以使用来自各种库的大量方法，例如LAPACK或BLAS，具体取决于矩阵A(正定,...)的属性。作为起点，我建议一种直接反演矩阵 A 的简单方法，尤其是在矩阵很小的情况下。还有许多迭代方法(雅可比，梯度，高斯赛德尔......)可用于更大的情况。

编辑：一个想法，分两步解决

第一步：您可以在 2 个子系统中重写您的系统，这些子系统将 X 和 Y 作为未知数，但维度等于每个向量中的未知数。 X 中的第一个子系统将是 AX = b，可以通过直接或迭代方法求解。

第二步：一旦知道X，Y中的第二个系统可以直接求解，因为Y将以Y = A'X + b'的形式表示

我认为这种方法更笼统。