我读过的关于该主题的最好的书是 Henk A. van der Vorst 的大型线性系统的迭代克雷洛夫方法。令人高兴的是，这本书的价格很合理。当然，在解决范德沃斯特问题之前，您需要在一般线性代数技术方面有坚实的基础。我推荐乔尔·N·富兰克林 (Joel N. Franklin) 1968 年出版的短书《矩阵理论》(Matrix Theory)，该书 1993 年以廉价平装本形式从多佛 (Dover) 购买。 （Van der Vorst 并没有提到 Kapp 和 Brown 的巧妙、简单且通常有效的有序多重交互方法，我多次发现该方法很有用，所以让我在这里提一下，以便您需要时可以查找。披露：我个人认识布朗，因此可能会偏袒；但是，他和卡普的技术仍然是一种很好的技术，范德沃斯特确实省略了。）
我和我的同事都尝试过 LAPACK 和 BLAS 的各种 C 和 C++ 接口。由于各种原因，我们发现他们中没有一个是完全令人满意的，尽管他们确实试图提供帮助。我们每个人最终决定完全跳过接口并直接使用 LAPACK 和 BLAS。这就是我倾向于向您推荐的。
LAPACK 和 BLAS 应该从 FORTRAN 77 中调用，而不是从 C++ 中调用。尽管如此，从 C++ 调用它们并不是那么困难，而且您不需要 C++ 接口来完成它。事实上，您可能不想要 C++ 接口，至少不想要其他人准备的通用接口。如果你不把 LAPACK 和 BLAS 合二为一的话，他们会更高兴。 （请记住，FORTRAN 77 的链接功能虽然有效，但缺乏 C 风格的头文件。）
要直接使用 LAPACK 和 BLAS，您必须知道的第一件事就是这一点，除非您理解它，否则您会很痛苦：矩阵按列优先存储。也就是说，矩阵的每一列（而不是每一行）在内存中顺序表示为一个单元。因此，矩阵元素直接存储在其上方和下方的元素之间，而不是直接存储在其左侧和右侧的元素之间。事实上，LAPACK和BLAS以这种方式存储矩阵是谨慎的，因为发明C的Kernighan和Ritchie似乎对线性代数从来没有很感兴趣。 C 是一种很好的语言，但 C 的默认矩阵存储约定可能从一开始就是一个错误。从数学家的角度来看，LAPACK 和 BLAS 以应有的存储方式存储矩阵，即列优先；而且，如果您正在编写线性代数代码，那么您也应该以这种方式存储矩阵。忽略 C 的默认行优先约定：它与矩阵数学的约定无关，并且您不需要它。
下面是一个完整的例子：
#include <iostream>
extern "C" {
void dgesv_(
const int *N,
const int *NRHS,
double *A,
const int *LDA,
int *IPIV,
double *B,
const int *LDB,
int *INFO
);
}
int main()
{
const int N = 2;
// Let A = | 3.0 6.1 |
// |-1.2 1.7 |
double A[N*N] = {3.0, -1.2, 6.1, 1.7};
// Let b = | 5.5 |
// | 0.4 |
double x[N] = {5.5, 0.4};
// (Why is b named x? Answer: because b and x share
// the same storage, because Lapack will write x over b.)
// Solve the linear system A*x = b.
const int NRHS = 1;
const int LDA = N;
int IPIV [N];
const int LDB = N;
int INFO;
dgesv_(&N, &NRHS, A, &LDA, IPIV, x, &LDB, &INFO);
// Uncomment the next line if you wish to see
// the error code Lapack returns.
//std::cout << "INFO == " << INFO << "\n";
// Output x.
for (int i = 0; i < N; ++i) std::cout << x[i] << "\n";
return 0;
}
显然，如果您正在寻找快速解决方案，我没有。可能有一个预装软件包可以完全满足您的要求，但我的经验表明，这种可能性对您不利。实践中会出现数百种不同类型的矩阵问题，预装软件无法破解。然而，LAPACK 和 BLAS 以及 Franklin、van der Vorst 以及您现在正在阅读的答案，提供了解决您的特定问题所需的所有工具。
祝你好运。
Though you may not realize it, your question is probably more a linear-algebra question than a programming question, though it does have a programming dimension to it.
In this answer, I assume that your actual problem could be much larger than the N = 8 sketch you have posed in illustration. The mathematical discussion comes first. Sample code follows.
The nature of linear-algebra problems seems to be that they are either practically treatable by general linear-algebra techniques or they are not. General linear-algebra techniques such as the L-U factorization and the singular-value decomposition have been thoroughly developed since the 1960s, and have a firm theoretical foundation. The trouble is that such general techniques tend to O(N*N*N) in computational treatment, which is to say that making your vector 10 times longer wants 1000 times the computer time to treat. Beyond some maximum practical size, your problem grows computationally intractable -- at least by general linear algebra techniques.
IF YOUR PROBLEM IS EXTREMELY LARGE
Fortunately, there exist a bewildering variety of special-purpose techniques to reduce the O(N*N*N) problem to an O(N*N) problem. These techniques have mostly emerged since the 1960s and are an area of active research. Unfortunately, to know which technique to apply, and how, requires an intimate understanding of the problem and a fairly good grasp of matrix mathematics. At O(N*N*N), general techniques are known. At O(N*N), there are no general techniques, only lots of special-purpose techniques that work on restricted kinds of systems. You can often find a suitable restriction, if you work at it long enough, but it is never as simple as dumping your matrix into some general-purpose solver. Not at O(N*N).
The best book I have read on the subject is Henk A. van der Vorst's Iterative Krylov Methods for Large Linear Systems. Happily, the book's price is reasonable. Of course, you will need a firm grounding in general linear-algebra techniques before tackling van der Voorst. I suggest Joel N. Franklin's short 1968 book Matrix Theory, available in a cheap 1993 paperback from Dover. (Van der Vorst does not happen to mention Kapp's and Brown's clever, simple and often effective Method of Ordered Multiple Interactions, which I have repeatedly found useful, so let me mention that here, so that you can look it up if you need it. Disclosure: I know Brown personally and could be partial for this reason; but, still, his and Kapp's technique is a fine one which van der Vorst does omit.)
For reasons I do not fully understand, most of the recent techniques (though not Kapp's and Brown's) seem to prefer to work in real numbers, treating complex numbers as a special case or ignoring them entirely. You have not mentioned whether your numbers are complex, but if they are, then this may limit your options somewhat.
As an intermediate level between Franklin and van der Vorst, if you lack the time or interest to tackle the latter, you should at least look up and acquaint yourself with the method of conjugate gradients.
IF YOUR PROBLEM IS ONLY MODERATELY LARGE
If your problem is only moderately large -- say, N < 20,000 or so -- then your task is much easier. Forget van der Vorst in this case: you don't need him. Depending on whether you want your answer in milliseconds, seconds, minutes or hours (you haven't mentioned which, but it only affects the practical limit on N), you can tolerate O(N*N*N) performance, and in this case Franklin's general techniques from the 1960s are quite robust.
Fast, efficient, thorough and correct, the LAPACK library, and its low-level helper BLAS, are the standard in this area. You should use them.
LAPACK, BLAS AND C++
A colleage of mine and I have each tried various C and C++ interfaces to LAPACK and BLAS. For various reasons, we have found none of them wholly satisfactory, though they do try to be helpful. Each of us ultimately decided to skip the interfaces altogether and to use LAPACK and BLAS directly. This is what I tend to recommend to you.
LAPACK and BLAS were meant to be called from FORTRAN 77, not from C++. Still, calling them from C++ is not so hard, and you don't need a C++ interface to do it. In fact, you probably don't want a C++ interface, at least not a generic one prepared by somebody else. LAPACK and BLAS will be happier if you don't smother them in one. (Remember that FORTRAN 77, whose linking facilities though effective were limited, lacked C-style header files.)
The first thing you must know to use LAPACK and BLAS directly is this, and you will be miserable until you understand it: matrices are stored column-major. That is, each column, not each row, of the matrix is represented in memory sequentially as a unit. Thus, a matrix element is stored immediately between the elements above and below it, not immediately between the elements left and right of it. Actually, LAPACK and BLAS are prudent to store matrices this way, for Kernighan and Ritchie, who invented C, never seemed very interested in linear algebra. C is a fine language, but C's default matrix storage convention was probably a mistake from the first. From the mathematician's perspective, LAPACK and BLAS store matrices the way they ought to be stored, column-major; and, if you are writing linear-algebra code then you should store your matrices this way, too. Ignore C's default, row-major convention: it has nothing to do with the conventions of matrix mathematics and you don't want it.
Here is a complete example:
#include <iostream>
extern "C" {
void dgesv_(
const int *N,
const int *NRHS,
double *A,
const int *LDA,
int *IPIV,
double *B,
const int *LDB,
int *INFO
);
}
int main()
{
const int N = 2;
// Let A = | 3.0 6.1 |
// |-1.2 1.7 |
double A[N*N] = {3.0, -1.2, 6.1, 1.7};
// Let b = | 5.5 |
// | 0.4 |
double x[N] = {5.5, 0.4};
// (Why is b named x? Answer: because b and x share
// the same storage, because Lapack will write x over b.)
// Solve the linear system A*x = b.
const int NRHS = 1;
const int LDA = N;
int IPIV [N];
const int LDB = N;
int INFO;
dgesv_(&N, &NRHS, A, &LDA, IPIV, x, &LDB, &INFO);
// Uncomment the next line if you wish to see
// the error code Lapack returns.
//std::cout << "INFO == " << INFO << "\n";
// Output x.
for (int i = 0; i < N; ++i) std::cout << x[i] << "\n";
return 0;
}
[Why is the function named dgesv_()? Answer: it isn't, at least not in FORTRAN 77, in which it is named DGESV. However, FORTRAN 77 is case-insensitive and, when we translate that to C's linking conventions, we get dgesv_(). At least, that's what my colleague and I get on every Linux, BSD or OSX machine on which we've tried it. On Debian or Ubuntu, the shell command "readelf --symbols /usr/lib/liblapack.so | grep -i DGESV" discovers this. On a Microsoft platform, the C-linked symbol may have some other name: you'll have to investigate if this is relevant to you.]
LAPACK and BLAS are very good at what they do. You should use them.
Whether you store your data in C-style arrays (as the example does) or in C++-style std::valarrays is up to you. LAPACK and BLAS don't care, so long as the storage is column-major.
THE MINIMIZATION OF COEFFICIENTS
You have not given enough information for me to tell exactly what you wish to do with your coefficients, but one suspects that what you may need is the solution to an N-by-(2*N) underdetermined system. Such runs beyond the scope of Franklin's book but, if you have already absorbed Franklin's material, then my own notes on the pseudoinverse may help you.
Obviously, if you were looking for a quick solution, I don't have it. There may be a precanned software package that does precisely what you want, but my experience suggests that the odds are against you in this. There are just too many hundreds of distinct kinds of matrix problems that arise in practice for precanned software to crack. However, LAPACK and BLAS, together with Franklin, van der Vorst and the answer you are now reading, offer all the tools you should need to crack your specific problem.
发布评论
评论（1）
尽管您可能没有意识到，您的问题可能更多是线性代数问题，而不是编程问题，尽管它确实具有编程维度。
在这个答案中，我假设您的实际问题可能比您在插图中提出的 N = 8 草图大得多。数学讨论是第一位的。示例代码如下。
线性代数问题的本质似乎是它们要么可以通过一般线性代数技术实际处理，要么不能。 LU 分解和奇异值分解等通用线性代数技术自 20 世纪 60 年代以来已得到彻底发展，并具有坚实的理论基础。问题在于，这种通用技术在计算处理上往往是 O(N*N*N)，也就是说，让你的向量长 10 倍需要 1000 倍的计算机处理时间。超出某个最大实际大小，您的问题在计算上变得难以解决——至少通过一般线性代数技术是这样。
如果您的问题非常大
幸运的是，存在多种令人眼花缭乱的专用技术可以将 O(N*N*N) 问题简化为 O(N*N) 问题。这些技术大多自 20 世纪 60 年代以来出现，是一个活跃的研究领域。不幸的是，要知道应用哪种技术以及如何应用，需要对问题有深入的了解并且对矩阵数学有相当好的掌握。在 O(N*N*N) 下，通用技术是已知的。在 O(N*N) 中，没有通用技术，只有许多适用于受限类型系统的专用技术。如果您工作足够长的时间，您通常可以找到合适的限制，但它永远不会像将矩阵转储到某些通用求解器中那么简单。不是 O(N*N)。
我读过的关于该主题的最好的书是 Henk A. van der Vorst 的大型线性系统的迭代克雷洛夫方法。令人高兴的是，这本书的价格很合理。当然，在解决范德沃斯特问题之前，您需要在一般线性代数技术方面有坚实的基础。我推荐乔尔·N·富兰克林 (Joel N. Franklin) 1968 年出版的短书《矩阵理论》(Matrix Theory)，该书 1993 年以廉价平装本形式从多佛 (Dover) 购买。 （Van der Vorst 并没有提到 Kapp 和 Brown 的巧妙、简单且通常有效的有序多重交互方法，我多次发现该方法很有用，所以让我在这里提一下，以便您需要时可以查找。披露：我个人认识布朗，因此可能会偏袒；但是，他和卡普的技术仍然是一种很好的技术，范德沃斯特确实省略了。）
出于我不完全理解的原因，大多数最近的技术（尽管不是卡普和布朗）似乎更喜欢使用实数，将复数视为特殊情况或完全忽略它们。您没有提到您的数字是否复杂，但如果是，那么这可能会在一定程度上限制您的选择。
作为富兰克林和范德沃斯特之间的中间水平，如果您缺乏时间或兴趣来解决后者，您至少应该查找并熟悉共轭梯度方法。
如果您的问题只是中等大
如果您的问题只是中等大 - 比如说，N << 20,000 左右——那么你的任务就容易多了。在这种情况下忘记范德沃斯特：你不需要他。根据你想要的答案是毫秒、秒、分钟还是小时（你没有提到哪一个，但它只影响 N 的实际限制），你可以容忍 O(N*N*N) 性能，在这个案例 富兰克林 20 世纪 60 年代的一般技术相当稳健。
LAPACK 库及其低级助手 BLAS 快速、高效、彻底且正确，是该领域的标准。你应该使用它们。
LAPACK、BLAS 和 C++
我和我的同事都尝试过 LAPACK 和 BLAS 的各种 C 和 C++ 接口。由于各种原因，我们发现他们中没有一个是完全令人满意的，尽管他们确实试图提供帮助。我们每个人最终决定完全跳过接口并直接使用 LAPACK 和 BLAS。这就是我倾向于向您推荐的。
LAPACK 和 BLAS 应该从 FORTRAN 77 中调用，而不是从 C++ 中调用。尽管如此，从 C++ 调用它们并不是那么困难，而且您不需要 C++ 接口来完成它。事实上，您可能不想要 C++ 接口，至少不想要其他人准备的通用接口。如果你不把 LAPACK 和 BLAS 合二为一的话，他们会更高兴。 （请记住，FORTRAN 77 的链接功能虽然有效，但缺乏 C 风格的头文件。）
要直接使用 LAPACK 和 BLAS，您必须知道的第一件事就是这一点，除非您理解它，否则您会很痛苦：矩阵按列优先存储。也就是说，矩阵的每一列（而不是每一行）在内存中顺序表示为一个单元。因此，矩阵元素直接存储在其上方和下方的元素之间，而不是直接存储在其左侧和右侧的元素之间。事实上，LAPACK和BLAS以这种方式存储矩阵是谨慎的，因为发明C的Kernighan和Ritchie似乎对线性代数从来没有很感兴趣。 C 是一种很好的语言，但 C 的默认矩阵存储约定可能从一开始就是一个错误。从数学家的角度来看，LAPACK 和 BLAS 以应有的存储方式存储矩阵，即列优先；而且，如果您正在编写线性代数代码，那么您也应该以这种方式存储矩阵。忽略 C 的默认行优先约定：它与矩阵数学的约定无关，并且您不需要它。
下面是一个完整的例子：
【为什么函数名为 dgesv_()？答案：不是，至少在 FORTRAN 77 中不是，它被命名为 DGESV。然而，FORTRAN 77 不区分大小写，当我们将其转换为 C 的链接约定时，我们得到 dgesv_()。至少，这是我和我的同事在我们尝试过的每台 Linux、BSD 或 OSX 机器上得到的结果。在 Debian 或 Ubuntu 上，shell 命令“readelf --symbols /usr/lib/liblapack.so | grep -i DGESV”会发现这一点。在 Microsoft 平台上，C 连接符号可能有其他名称：您必须调查这是否与您相关。]
LAPACK 和 BLAS 非常擅长它们的工作。你应该使用它们。
将数据存储在 C 样式数组（如示例所示）中还是 C++ 样式 std::valarrays 中取决于您。 LAPACK 和 BLAS 并不关心，只要存储是列优先的即可。
系数的最小化
您没有提供足够的信息让我准确地告诉我您希望如何处理系数，但有人怀疑您可能需要的是 N×(2*N) 欠定系统的解。这超出了富兰克林书的范围，但是，如果您已经吸收了富兰克林的材料，那么我自己关于伪逆的笔记可能会帮助你。
显然，如果您正在寻找快速解决方案，我没有。可能有一个预装软件包可以完全满足您的要求，但我的经验表明，这种可能性对您不利。实践中会出现数百种不同类型的矩阵问题，预装软件无法破解。然而，LAPACK 和 BLAS 以及 Franklin、van der Vorst 以及您现在正在阅读的答案，提供了解决您的特定问题所需的所有工具。
祝你好运。
Though you may not realize it, your question is probably more a linear-algebra question than a programming question, though it does have a programming dimension to it.
In this answer, I assume that your actual problem could be much larger than the N = 8 sketch you have posed in illustration. The mathematical discussion comes first. Sample code follows.
The nature of linear-algebra problems seems to be that they are either practically treatable by general linear-algebra techniques or they are not. General linear-algebra techniques such as the L-U factorization and the singular-value decomposition have been thoroughly developed since the 1960s, and have a firm theoretical foundation. The trouble is that such general techniques tend to O(N*N*N) in computational treatment, which is to say that making your vector 10 times longer wants 1000 times the computer time to treat. Beyond some maximum practical size, your problem grows computationally intractable -- at least by general linear algebra techniques.
IF YOUR PROBLEM IS EXTREMELY LARGE
Fortunately, there exist a bewildering variety of special-purpose techniques to reduce the O(N*N*N) problem to an O(N*N) problem. These techniques have mostly emerged since the 1960s and are an area of active research. Unfortunately, to know which technique to apply, and how, requires an intimate understanding of the problem and a fairly good grasp of matrix mathematics. At O(N*N*N), general techniques are known. At O(N*N), there are no general techniques, only lots of special-purpose techniques that work on restricted kinds of systems. You can often find a suitable restriction, if you work at it long enough, but it is never as simple as dumping your matrix into some general-purpose solver. Not at O(N*N).
The best book I have read on the subject is Henk A. van der Vorst's Iterative Krylov Methods for Large Linear Systems. Happily, the book's price is reasonable. Of course, you will need a firm grounding in general linear-algebra techniques before tackling van der Voorst. I suggest Joel N. Franklin's short 1968 book Matrix Theory, available in a cheap 1993 paperback from Dover. (Van der Vorst does not happen to mention Kapp's and Brown's clever, simple and often effective Method of Ordered Multiple Interactions, which I have repeatedly found useful, so let me mention that here, so that you can look it up if you need it. Disclosure: I know Brown personally and could be partial for this reason; but, still, his and Kapp's technique is a fine one which van der Vorst does omit.)
For reasons I do not fully understand, most of the recent techniques (though not Kapp's and Brown's) seem to prefer to work in real numbers, treating complex numbers as a special case or ignoring them entirely. You have not mentioned whether your numbers are complex, but if they are, then this may limit your options somewhat.
As an intermediate level between Franklin and van der Vorst, if you lack the time or interest to tackle the latter, you should at least look up and acquaint yourself with the method of conjugate gradients.
IF YOUR PROBLEM IS ONLY MODERATELY LARGE
If your problem is only moderately large -- say, N < 20,000 or so -- then your task is much easier. Forget van der Vorst in this case: you don't need him. Depending on whether you want your answer in milliseconds, seconds, minutes or hours (you haven't mentioned which, but it only affects the practical limit on N), you can tolerate O(N*N*N) performance, and in this case Franklin's general techniques from the 1960s are quite robust.
Fast, efficient, thorough and correct, the LAPACK library, and its low-level helper BLAS, are the standard in this area. You should use them.
LAPACK, BLAS AND C++
A colleage of mine and I have each tried various C and C++ interfaces to LAPACK and BLAS. For various reasons, we have found none of them wholly satisfactory, though they do try to be helpful. Each of us ultimately decided to skip the interfaces altogether and to use LAPACK and BLAS directly. This is what I tend to recommend to you.
LAPACK and BLAS were meant to be called from FORTRAN 77, not from C++. Still, calling them from C++ is not so hard, and you don't need a C++ interface to do it. In fact, you probably don't want a C++ interface, at least not a generic one prepared by somebody else. LAPACK and BLAS will be happier if you don't smother them in one. (Remember that FORTRAN 77, whose linking facilities though effective were limited, lacked C-style header files.)
The first thing you must know to use LAPACK and BLAS directly is this, and you will be miserable until you understand it: matrices are stored column-major. That is, each column, not each row, of the matrix is represented in memory sequentially as a unit. Thus, a matrix element is stored immediately between the elements above and below it, not immediately between the elements left and right of it. Actually, LAPACK and BLAS are prudent to store matrices this way, for Kernighan and Ritchie, who invented C, never seemed very interested in linear algebra. C is a fine language, but C's default matrix storage convention was probably a mistake from the first. From the mathematician's perspective, LAPACK and BLAS store matrices the way they ought to be stored, column-major; and, if you are writing linear-algebra code then you should store your matrices this way, too. Ignore C's default, row-major convention: it has nothing to do with the conventions of matrix mathematics and you don't want it.
Here is a complete example:
[Why is the function named dgesv_()? Answer: it isn't, at least not in FORTRAN 77, in which it is named DGESV. However, FORTRAN 77 is case-insensitive and, when we translate that to C's linking conventions, we get dgesv_(). At least, that's what my colleague and I get on every Linux, BSD or OSX machine on which we've tried it. On Debian or Ubuntu, the shell command "readelf --symbols /usr/lib/liblapack.so | grep -i DGESV" discovers this. On a Microsoft platform, the C-linked symbol may have some other name: you'll have to investigate if this is relevant to you.]
LAPACK and BLAS are very good at what they do. You should use them.
Whether you store your data in C-style arrays (as the example does) or in C++-style std::valarrays is up to you. LAPACK and BLAS don't care, so long as the storage is column-major.
THE MINIMIZATION OF COEFFICIENTS
You have not given enough information for me to tell exactly what you wish to do with your coefficients, but one suspects that what you may need is the solution to an N-by-(2*N) underdetermined system. Such runs beyond the scope of Franklin's book but, if you have already absorbed Franklin's material, then my own notes on the pseudoinverse may help you.
Obviously, if you were looking for a quick solution, I don't have it. There may be a precanned software package that does precisely what you want, but my experience suggests that the odds are against you in this. There are just too many hundreds of distinct kinds of matrix problems that arise in practice for precanned software to crack. However, LAPACK and BLAS, together with Franklin, van der Vorst and the answer you are now reading, offer all the tools you should need to crack your specific problem.
Good luck.