jjzjj

c++ - 非常小方阵的特征线性求解器

coder 2024-02-04 原文

我正在使用 Eigen一个求解极小方阵(4X4)线性方程的C++程序。

我的测试代码是这样的

template<template <typename MatrixType> typename EigenSolver>
Vertor3d solve(){
   //Solve Ax = b and A is a real symmetric matrix and positive semidefinite

   ... // Construct 4X4 square matrix A and 4X1 vector b

   EigenSolver<Matrix4d> solver(A);
   auto x = solver.solve(b);

   ... // Compute relative error for validating
}

我测试了一些 EigenSolver其中包括:

  1. FullPixLU
  2. PartialPivLU
  3. HouseholderQR
  4. ColPivHouseholderQR
  5. ColPivHouseholderQR
  6. 完全正交分解
  7. 低密度脂蛋白
  8. 正反

正逆是:

template<typename MatrixType>
struct InverseSolve
{
private:
    MatrixType  inv;
public:
    InverseSolve(const MatrixType &matrix) :inv(matrix.inverse()) {
    }
    template<typename VectorType>
    auto solve(const VectorType & b) {
        return inv * b;
    }
};

我发现最快的方法是DirectInverse，即使我把Eigen和MKL连接起来，结果也没有改变。

这是测试结果

  1. FullPixLU:477 毫秒
  2. PartialPivLU:468 毫秒
  3. HouseholderQR:849 毫秒
  4. ColPivHouseholderQR:766 毫秒
  5. ColPivHouseholderQR:857 毫秒
  6. 完全正交分解:832 毫秒
  7. LDLT:477 毫秒
  8. 正反:88 毫秒

它全部使用 1000000 个矩阵，从均匀分布 [0,100] 随机加倍。我首先构造上三角，然后复制到下三角。

DirectInverse 唯一的问题是它的相对误差比其他求解器略大但可以接受。

我的程序是否有更快或更优雅的解决方案？DirectInverse 是我程序的快速解决方案吗？

DirectInverse 不使用对称信息，为什么 DirectInverse 比 LDLT 快得多？

最佳答案

尽管许多人建议当您只想求解线性系统时永远不要显式计算逆，但对于非常小的矩阵，这实际上是有益的，因为存在使用余因子的封闭形式的解决方案。

您测试的所有其他替代方案都会变慢，因为它们会进行旋转(这意味着分支)，即使对于小型固定大小的矩阵也是如此。此外，它们中的大多数会导致更多的划分，并且不像直接计算那样可向量化。

为了提高准确性(如果需要，该技术实际上可以独立于求解器使用)，您可以通过使用残差再次求解系统来改进初始解:

Eigen::Vector4d solveDirect(const Eigen::Matrix4d& A, const Eigen::Vector4d& b)
{
    Eigen::Matrix4d inv = A.inverse();
    Eigen::Vector4d x = inv * b;
    x += inv*(b-A*x);
    return x;
}

我不认为 Eigen 在这里直接提供了一种利用 A 对称性的方法(对于直接计算的逆)。您可以尝试通过将 A 的自伴 View 显式复制到临时文件中来暗示这一点，并希望编译器足够聪明以找到公共(public)子表达式:

Eigen::Matrix4d tmp = A.selfadjointView<Eigen::Upper>();
Eigen::Matrix4d inv = tmp.inverse();

为了减少一些除法，您还可以使用 -freciprocal-math(在 gcc 或 clang 上)进行编译，这当然会稍微降低准确性。

如果这确实对性能至关重要，请尝试实现手动调整的 inverse_4x4_symmetric 方法。 利用 inv * b 的对称性不太可能对如此小的矩阵有益。

关于c++ - 非常小方阵的特征线性求解器，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/50909385/

小方求解codeli毫秒c++linear-algebraeigen

有关c++ - 非常小方阵的特征线性求解器的更多相关文章

  1. ruby-on-rails - 如何优雅地重启 thin + nginx？ - 2

    我的瘦服务器配置了nginx，我的ROR应用程序正在它们上运行。在我发布代码更新时运行thinrestart会给我的应用程序带来一些停机时间。我试图弄清楚如何优雅地重启正在运行的Thin实例，但找不到好的解决方案。有没有人能做到这一点？ 最佳答案 #Restartjustthethinserverdescribedbythatconfigsudothin-C/etc/thin/mysite.ymlrestartNginx将继续运行并代理请求。如果您将Nginx设置为使用多个上游服务器，例如server{listen80;server

  2. ruby - 如何根据特征实现 FactoryGirl 的条件行为 - 2

    我有一个用户工厂。我希望默认情况下确认用户。但是鉴于unconfirmed特征，我不希望它们被确认。虽然我有一个基于实现细节而不是抽象的工作实现，但我想知道如何正确地做到这一点。factory:userdoafter(:create)do|user,evaluator|#unwantedimplementationdetailshereunlessFactoryGirl.factories[:user].defined_traits.map(&:name).include?(:unconfirmed)user.confirm!endendtrait:unconfirmeddoenden

  3. ruby - 使用 `+=` 和 `send` 方法 - 2

    如何将send与+=一起使用？a=20;a.send"+=",10undefinedmethod`+='for20:Fixnuma=20;a+=10=>30 最佳答案 恐怕你不能。+=不是方法，而是语法糖。参见http://www.ruby-doc.org/docs/ProgrammingRuby/html/tut_expressions.html它说Incommonwithmanyotherlanguages,Rubyhasasyntacticshortcut:a=a+2maybewrittenasa+=2.你能做的最好的事情是:

  4. ruby-on-rails - 使用 HTTParty 的非常基本的 Rails 4.1 API 调用 - 2

    Rails相对较新。我正在尝试调用一个API，它应该向我返回一个唯一的URL。我的应用程序中捆绑了HTTParty。我已经创建了一个UniqueNumberController，并且我已经阅读了几个HTTParty指南，直到我想要什么，但也许我只是有点迷路，真的不知道该怎么做。基本上，我需要做的就是调用API，获取它返回的URL，然后将该URL插入到用户的数据库中。谁能给我指出正确的方向或与我分享一些代码？ 最佳答案 假设API为JSON格式并返回如下数据:{"url":"http://example.com/unique-url"

  5. ruby - 如何计算 Liquid 中的变量 +1 - 2

    我对如何计算通过{%assignvar=0%}赋值的变量加一完全感到困惑。这应该是最简单的任务。到目前为止，这是我尝试过的:{%assignamount=0%}{%forvariantinproduct.variants%}{%assignamount=amount+1%}{%endfor%}Amount:{{amount}}结果总是0。也许我忽略了一些明显的东西。也许有更好的方法。我想要存档的只是获取运行的迭代次数。 最佳答案 因为{{incrementamount}}将输出您的变量值并且不会影响{%assign%}定义的变量，我

  6. ruby - cucumber 特征和步骤定义 - 2

    我是Cucumber测试的新手。我创建了两个特征文件:events.featurepartner.feature并将我的步骤定义放在step_definitions文件夹中:./step_definitions/events.rbpartner.rbCucumber似乎在所有.rb文件中查找步骤信息。有没有办法限制该功能查看特定的步骤定义文件？我之所以要这样做，是因为即使我使用了--guess标志，我也会遇到不明确的匹配错误。我之所以要这样做，有以下几个原因。我正在测试CMS，并希望在不同的功能中测试每种不同的内容类型(事件和合作伙伴)。事件.特征Feature:AddpartnerA

  7. arrays - Ruby 数组 += vs 推送 - 2

    我有一个数组数组，想将元素附加到子数组。+=做我想做的，但我想了解为什么push不做。我期望的行为(并与+=一起工作):b=Array.new(3,[])b[0]+=["apple"]b[1]+=["orange"]b[2]+=["frog"]b=>[["苹果"],["橙子"],["Frog"]]通过推送，我将推送的元素附加到每个子数组(为什么？):a=Array.new(3,[])a[0].push("apple")a[1].push("orange")a[2].push("frog")a=>[[“苹果”、“橙子”、“Frog”]、[“苹果”、“橙子”、“Frog”]、[“苹果”、“

  8. += 的 Ruby 方法 - 2

    有没有办法让Ruby能够做这样的事情？classPlane@moved=0@x=0defx+=(v)#thisiserror@x+=v@moved+=1enddefto_s"moved#{@moved}times,currentxis#{@x}"endendplane=Plane.newplane.x+=5plane.x+=10putsplane.to_s#moved2times,currentxis15 最佳答案 您不能在Ruby中覆盖复合赋值运算符。任务在内部处理。您应该覆盖+，而不是+=。plane.a+=b与plane.a=

  9. ruby - 如何在 Ruby 中生成一个非常大的随机整数？ - 2

    我想在ruby​​中生成一个64位整数。我知道在Java中你有很多渴望，但我不确定你会如何在Ruby中做到这一点。另外，64位数字中有多少个字符？这是我正在谈论的示例......123456789999。@num=Random.rand(9000)+Random.rand(9000)+Random.rand(9000)但我认为这是非常低效的，必须有一种更简单、更简洁的方法来做到这一点。谢谢! 最佳答案 rand可以将范围作为参数:pa=rand(2**32..2**64-1)#=>11093913376345012184putsa.

  10. ruby - Sinatra + Heroku + Datamapper 使用 dm-sqlite-adapter 部署问题 - 2

    出于某种原因，heroku尝试要求dm-sqlite-adapter，即使它应该在这里使用Postgres。请注意，这发生在我打开任何URL时-而不是在gitpush本身期间。我构建了一个默认的Facebook应用程序。gem文件:source:gemcuttergem"foreman"gem"sinatra"gem"mogli"gem"json"gem"httparty"gem"thin"gem"data_mapper"gem"heroku"group:productiondogem"pg"gem"dm-postgres-adapter"endgroup:development,:t

随机推荐