拓冰网站定制 
简介面向 Visual C 开发者的 newmat 矩阵运算库压缩包能在 C 项目中快速完成矩阵加减乘、转置、求逆、特征值等数值计算适合科学计算、图像处理与信号处理等需要密集矩阵运算的场景。包内共 97 个文件大小约 284KB以 56 个 cpp 源文件、12 个 h 头文件为主体另含工程文件、txt 说明文档和示例程序既有库实现也有可运行的测试例程便于对照学习与二次封装。目前已有 163 人学习使用。对希望摆脱底层矩阵算法细节、专注业务建模的 C 程序员来说这是一套开箱即用的数学计算工具能有效提升项目开发效率。1. newmat.lib 放进 VC 工程里不只是“能算”拿到这份 newmat.lib.zip 时第一反应是“又是一个被打包好的旧库”。但把压缩包翻开后才注意到里面不仅有newmat.lib这个编译产物和newmat.h主头文件还包括newmat1.cpp到newmat9.cpp的全部实现源码以及cholesky.cpp、svd.cpp、jacobi.cpp、hholder.cpp、evalue.cpp、newfft.cpp这些专门做数值算法的小模块。对用 Visual C 下手的人来说这意味着两件事一是不用为 BLAS/LAPACK 的环境挨个配置二是不必依赖别人打包那一个.lib的编译环境源码在手随时能重编出匹配自己工程 ABI 的静态库。接下来会把库的源码结构、VC 集成方式、算法模块边界和打包示例代码都过一遍每个步骤说明为什么这么设。2. 拆包后先分清源码构成再决定用哪条路径进 VS 工程2.1 三类文件别混编先建立映射解压后文件很多先按角色划分。newmat.h是主头文件newmatap.h是高级数值接口newmatio.h提供输出include.h在编译期调整 bool、STL 等兼容性设置。实现层分两组newmat1.cpp到newmat9.cpp加上bandmat.cpp、submat.cpp、myexcept.cpp是基础矩阵运算cholesky.cpp、svd.cpp、jacobi.cpp、hholder.cpp、evalue.cpp、newfft.cpp、sort.cpp、nm_misc.cpp、newmatnl.cpp是数值算法。剩下tmt*.cpp、nm_ex*.cpp、garch.cpp、nl_ex.cpp、sl_ex.cpp是测试和示例不要混进库工程。角色主要文件在 Visual Studio 工程里的放置基础类型头文件include.h, newmat.h, newmatrc.h, myexcept.hinclude 目录高级算法头文件newmatap.h, newmatnl.hinclude 目录核心实现newmat1.cpp~newmat9.cpp, bandmat.cpp, submat.cpp, myexcept.cpp编入静态库数值算法实现cholesky.cpp, svd.cpp, jacobi.cpp, hholder.cpp, evalue.cpp, newfft.cpp, nm_misc.cpp, newmatnl.cpp按需编入静态库测试与示例tmt1.cpp~tmtm.cpp, nm_ex1/2/3.cpp, garch.cpp, nl_ex.cpp, sl_ex.cpp独立可执行工程不进库这份映射是我在多个用 newmat 的老工程里验证过的。最容易踩的问题是把tmt1.cpp也加进库工程然后发现符号重复和额外的 main 入口。库工程只产出静态库测试代码单独生成 exe。2.2 Visual Studio 工程方式编出供自己用的 newmat.lib在 VS 里新建“静态库”项目目标平台选 x64。把上一节第二、三行所有实现文件拖进“源文件”头文件放“头文件”目录。配置 Debug/Release 两套都编译要特别注意语言选项设置“C/C ▸ 代码生成”里的运行时库为“多线程 DLL (/MD)”这样使用方工程无论走 Release 还是 Debug 都能链上对应的 CRT。异常模型保持/EHscnewmat 内部大量依赖异常上报越界。编译选项做一句话定位/O2开启速度优化newmat 这类数值库在 O2 下运行效率接近普通手写裸循环/MD是大多数现代应用项目的选择如果换成/MT静态 CRT 也不是不能跑但后续链接另一个也用了/MT的第三方库时会频繁报 LNK4098 警告。如果最终工程使用 VS2015 之后的版本建议把库工程和使用方工程的“平台工具集”保持在同一版本避免运行时行为不一致。2.3 命令行快速重编一次CI 里也能用没有打开 IDE 时包里的nm_i8.mak、nm_b55.mak、nm_gnu.mak、nm_ow.mak分别对应 Intel、Borland、GCC、OpenWatcom但没提供一个现成的 MSVC makefile。我一般在 CI 脚本里直接调cl和lib两条命令cl /nologo /O2 /EHsc /MD /c newmat1.cpp newmat2.cpp newmat3.cpp newmat4.cpp \ newmat5.cpp newmat6.cpp newmat7.cpp newmat8.cpp newmat9.cpp bandmat.cpp \ submat.cpp myexcept.cpp cholesky.cpp svd.cpp jacobi.cpp hholder.cpp \ evalue.cpp newfft.cpp sort.cpp nm_misc.cpp newmatnl.cpp lib /NOLOGO /OUT:newmat.lib *.obj第一行把矩阵库全部实现源码分别编译成对象文件。参数含义很直接/c代表只编译不链接免得cl因为入口函数缺失报 LNK2019/O2对应速度优化对大型矩阵乘法收益较大/MD使用动态 CRT和主工程保持一致。第二行的lib是 MSVC 自带的库管理器把上一行产出的所有.obj包进newmat.lib。整个过程不用配置 VS 工程文件之后无论是 CMake 还是 msbuild 项目都可以直接把newmat.lib当作外部依赖链接。提示如果工程同时链接别人给的一份 newmat.lib又把newmat1.cpp等源文件作为编译单元加进去链接阶段会出现复数符号定义错误。源码和 .lib 只能二选一。3. newmat 的矩阵类型与引用语义决定你该怎么写运算3.1 矩阵类型家族Matrix仍是中心但其余类型不是它的子类newmat 这里比较反直觉的一点是SymmetricMatrix、BandMat、DiagonalMatrix不是从Matrix继承的而是都从BaseMatrix这个抽象基类派生。这样做的目的是让重载运算符的返回值类型能按具体运算保持比如SymmetricMatrix SymmetricMatrix可以直接返回SymmetricMatrix如果强迫先转成Matrix再运算带宽和对称信息就浪费了。BandMat内部只存带内的lowerupper1行元素访问(i,j)时先检查是否落在带内越界时bandmat.cpp里的函数会抬一个 index error。写通用算法时参数要写成const BaseMatrix实际计算在Matrix构造时完成展开。常见的坑是直接拿const BaseMatrix调.t()、.i()看起来能用实际返回的是临时对象再用它参与运算容易把类型退化掉。所以我在团队里约定函数入参用const BaseMatrix但对入参做的最终运算必须在函数体内显式转成具体类型。运算对象类型存储内容常见误用Matrix全矩形数据行主序拿 0 当行下标SymmetricMatrix下三角存储不包括上三角写S(i,j)且i j之外的元素BandMat只存带宽内元素对带外元素赋值DiagonalMatrix单一主对角数组访问D(i,j)时i ! j3.2 运算符的临时对象机制不是所有表达式都能随意写newmat 内部做了完整的表达式临时对象复用A B C D这类长表达式能编译也能跑但如果你把同一矩阵放在赋值号两边例如A A B它能处理。真正危险的是A A A.t()这种在临时对象构造过程中复用了同一块存储的情况求转置时可能读到已经覆盖过的数据。安全示例是这样的Matrix A(2, 2); A 1.0 2.0 3.0 4.0; Matrix B A A.t(); // 安全临时对象类型明确 Matrix C; C A A.t(); // 也安全 Matrix D A * A; // 安全结果矩阵是独立分配的 label:这里B是拷贝构造C先默认构造再走赋值运算符两种情况都不共享数据。需要真正避免的是 in-place 的更新newmat 没有实现类似A A.t()这种自反语句优化对于SymmetricMatrix尤其容易出问题。安全做法是先用临时对象算出结果再赋回去。代码块的注释不用理会label:只是给断点用的一个空标签编译型代码里可以删掉。3.3 子矩阵获取columns() 返回的不是裸引用submat.cpp里实现了一套较完整的子矩阵接口A.rows(a, b)、A.columns(c, d)、A.submatrix(r1, r2, c1, c2)都返回可分派给Matrix的结果对象。这里很多人被“视图”两个字误导以为拿到的是不需要拷贝的引用。其实 newmat 的做法是返回一个持有原矩阵引用的临时对象如果把它直接赋给一个新的Matrixnewmat 会做一次完整复制。如果试图用非 const 引用去绑定这个临时对象编译器直接报错Matrix A(4, 4); Matrix M A.columns(1, 2); // 安全M 获得复制后的独立数据 // Matrix M2 A.columns(1, 2); // 错误临时对象不能绑定非 const 引用推荐做法是子矩阵赋给Matrix时就复制后续对M的任何修改不影响A。如果确实要原地更新一行用MatrixRow系列接口比如MatrixRow mr(A, 1); mr A.Row(2);这会把第二行的数据复制到第一行而不会产生新的矩阵对象。这样写语义清楚也符合 newmat 内部引用的生命周期。4. Cholesky、SVD、特征值与 FFTnewmatap 模块的调用边界4.1 Cholesky 分解构造、判断和反代求 solvecholesky.cpp对应newmatap.h里的Cholesky类。标准用法是构造时就完成分解。构造函数有两个参数第二个是可选的最小主元容忍度 epsilon一般保持默认即可。#include newmat.h #include newmatap.h SymmetricMatrix A(3); A 4.0 -2.0 -2.0 5.0 1.0 6.0; Cholesky chol(A); if (!chol.IsValid()) { cerr matrix not SPD endl; return -1; } ColumnVector b(3); b 8.0 0.0 7.0; ColumnVector x chol.i(b);在调用顺序上要提醒chol.i(b)不是先求逆矩阵再乘向量而是内部直接做 forward/back substitution得到线性方程组的解。随后取chol.cholesky()也能正常工作因为分解结果已缓存在对象内部。IsValid()的返回值基于分解过程中对角线元素与 epsilon 的比较判断所以一个理论上正定但条件数很差的矩阵可能被判无效这一点和 MATLAB 的 chol 默认行为不一样遇到边界矩阵时要有心理准备。4.2 SVD 与伪逆奇异值截断要自己处理svd.cpp提供SVD(Matrix, DiagonalMatrix, Matrix, Matrix)。它输出的是 V不是 VT。写最小二乘时我通常手动构造逆奇异值矩阵而不是调用D.i()Matrix U, V; DiagonalMatrix D; SVD(A, D, U, V); Real tol 1e-10; DiagonalMatrix Dinv(D.nrows()); for (int i 1; i D.nrows(); i) Dinv(i, i) (fabs(D(i, i)) tol) ? 1.0 / D(i, i) : 0.0; Matrix pseudoInv V * Dinv * U.t();这里Dinv(i, i)对小于 tol 的奇异值直接置 0而不是调D.i()否则 0 奇异值的倒数会变成一个异常大的数让最后结果退化。D.nrows()返回奇异值个数矩阵 U 和 V 的维度与此匹配。伪逆矩阵算出来后我再做一次pseudoInv * A的乘积检查确认结果接近单位阵避免在条件数很大的场景下自己没改干净。4.3 特征值Jacobi、Householder、EigenValues 的选用evalue.cpp、jacobi.cpp、hholder.cpp三个文件对应三套算法。实际使用中EigenValues是对称矩阵全特征值的首选接口SymmetricMatrix S(3); S 2.0 -1.0 0.0 2.0 -1.0 1.0; DiagonalMatrix D; Matrix V; EigenValues(S, D, V);D 中特征值按从大到小排序V 的每一列对应一个特征向量。如果数据量不大但需要高性能可用Jacobi做迭代如果矩阵维度较大且只求部分特征值先Householder化成三对角矩阵再对目标区间做二分和反迭代内存占用会低不少。这个选择顺序不是玄学Householder 的第一步是固定的正交变换代价是 O(n³)但之后求部分特征值只需要 O(n²) 量级的迭代整体比重启全量特征值快很多。4.4 FFT长度必须是 2 的幂否则别怪结果是 NaNnewfft.cpp里封装了FFT_Engine类。构造参数就是变换长度 n内部会按 n 构建旋转因子表。示例调用如下int n 8; FFT_Engine fft(n); ColumnVector signal(n), spectrum(n); for (int i 1; i n; i) signal(i) sin(2.0 * 3.14159265 * i / n); spectrum fft.Fft(signal);Fft返回频谱数据接口本身不复杂。复杂的是 n 的限制如果不是 2 的幂newfft.cpp内部的拆分逻辑会进入未定义分支轻则结果 NaN重则越界写内存。在调用前加一层检查把 n 向上对齐到不小于它的最小 2 的幂这是最稳妥的做法。任务类/函数限制求解对称正定线性方程组Cholesky必须对称正定最小二乘SVD需自行截断奇异值对称矩阵全部特征值EigenValues适用于对称矩阵对称矩阵三对角化Householder可用于部分特征值场景频谱分析FFT_Engine长度必须是 2 的幂5. 从 garch.cpp、nl_ex.cpp 到自己的统计模型示例代码的改造5.1 garch.cpp 是一个把矩阵求逆和循环放在一起的范例garch.cpp对多数人来说第一眼像统计模型但实际价值是演示怎样在既有项目里用 newmat 反复处理协方差矩阵和二次型。把它挪进 VC 工程照抄会先碰到两件事文件里大量用输出矩阵必须保证#include newmatio.h出现在所有输出之前另一个是文件使用Real类型在include.h里定义如果源码里直接写了double后续要全局替换。我一般把garch.cpp处理成目标函数只看它更新统计算子的部分#include newmat.h #include newmatap.h Real garch_log_likelihood(const ColumnVector y, const ColumnVector theta) { int T y.nrows(); Real omega theta(1); Real alpha theta(2); Real beta theta(3); Real h omega / (1.0 - alpha - beta); Real ll 0.0; for (int t 2; t T; t) { h omega alpha * y(t - 1) * y(t - 1) beta * h; ll -0.5 * (log(2.0 * 3.141592653589793) log(h) y(t) * y(t) / h); } return ll; }这里y是从garch.dat里读出的收益序列theta是三个待估参数。函数体本身没直接用矩阵运算但它验证了整个依赖链多变量 GARCH 的状态方程必然出现矩阵乘法和 Cholesky 求逆而这正是 newmat 的长处。zip 里同时放着nm_i10.mak、nm_m6.mak、nm_cc.mak也是提醒你这个库既可以被拆成单文件参与编译也可以交给多种 C 编译器。5.2 nl_ex.cpp 和 sl_ex.cpp不是玩具是调试工具链nl_ex.cpp演示的是NonLinearLeastSquares类把目标函数写成ReturnMatrix形式由库自己处理数值梯度和步长。这类代码看着是示例实际价值在于它提供了全局函数的调用范式当你在本地把nl_ex.cpp编译通过newmat 的异常捕获、临时矩阵管理、以及include.h里的编译器兼容开关就都已验证过。sl_ex.cpp则配合nm_misc.cpp里的Simplex函数给低维参数优化兜底。遇到目标函数梯度不稳定时先拿单纯形法跑一遍找初始值再用带梯度的优化器精修往往比直接上全局搜索收敛得更快。5.3 把示例改成可复用的封装只暴露业务接口从示例提炼规则时第一原则是不要在业务代码里散落 newmat 类型。一个最小封装是把 GARCH 的预测过程包成类newmat 类型全部留在 .cpp 内部class GarchModel { public: void setData(const std::vectordouble data); std::vectordouble predict(unsigned int horizon) const; private: ColumnVector r_; Real omega_, alpha_, beta_; }; void GarchModel::setData(const std::vectordouble data) { r_.ReSize(static_castint(data.size())); for (int i 0; i static_castint(data.size()); i) r_(i 1) data[i]; }注意r_.ReSize(n)不保证保留旧值与 STL 的resize行为不同。如果调用ReSize后想保留前几行旧数据需要先把旧值复制到临时数组。这个细节在实际更换数据源时踩过我调完ReSize后想保留前三行旧值结果全部归零。内部实现里predict再调用 5.1 的garch_log_likelihood配合优化器迭代外部业务只看到std::vector将来换掉底层库也不动接口。6. 用 tmt*.cpp 和 test_exc.cpp 验收库再谈三个旧库特有的坑6.1 先跑整套测试别急着量性能newmat 自带的tmt1.cpp到tmtm.cpp是针对不同模块的回归测试test_exc.cpp验证异常分支。把它们从库工程拆出来单独建一个 console 项目main 写在某个 tmt 文件里一次只编一个 tmt 文件逐个运行。最简单的验收命令是在 VS 开发命令行里做cl /nologo /O2 /EHsc newmat1.cpp newmat2.cpp newmat3.cpp newmat4.cpp \ newmat5.cpp newmat6.cpp newmat7.cpp newmat8.cpp newmat9.cpp bandmat.cpp \ submat.cpp myexcept.cpp cholesky.cpp svd.cpp jacobi.cpp hholder.cpp \ evalue.cpp newfft.cpp sort.cpp nm_misc.cpp newmatnl.cpp tmt1.cpp /Fe:tmt1.exe tmt1.exe正常输出末尾会出现类似 Ok 或 Tests completed 的信息。如果看到 fail多半不是库本身问题而是编译器对模板表达式的支持不同优先检查include.h顶部的条件编译宏。tmt 系列跑通后再有新的矩阵运算改动就能拿它们当基线回归集。6.2 三个旧库特有的坑绕不过去第一内存泄漏报告是性能问题不是真错误。newmat 自带 new/delete 统计控制台关闭前输出 Return from main - check for possible leaks很多初次使用者当失败处理。实际上只要连续两次运行泄漏数量不增长就可以认为正常。第二预编译的 newmat.lib 版本与主工程 CRT 不一致时异常流向会乱。压缩包里那个 newmat.lib 可能是多年前编译器编的和你当前 VS 的 C ABI 不匹配不要迷信它能直接链。拿到源码的价值就在于此重新执行第 2 节的cl和lib两条命令确保本地产物与主工程同一套工具链编译。第三USE_LAPACK宏。如果用源码编译却预定义了USE_LAPACK编译器会去找blas.h和lapack.h而包里并没有这两个头文件。检查include.h或工程预处理定义把没用的USE_LAPACK删掉很多 CI 报错都是因为这个宏从全局配置里传了下来。本文还有配套的精品资源点击获取 newmat库在Visual C++中的集成与数值算法实战指南
简介面向 Visual C 开发者的 newmat 矩阵运算库压缩包能在 C 项目中快速完成矩阵加减乘、转置、求逆、特征值等数值计算适合科学计算、图像处理与信号处理等需要密集矩阵运算的场景。包内共 97 个文件大小约 284KB以 56 个 cpp 源文件、12 个 h 头文件为主体另含工程文件、txt 说明文档和示例程序既有库实现也有可运行的测试例程便于对照学习与二次封装。目前已有 163 人学习使用。对希望摆脱底层矩阵算法细节、专注业务建模的 C 程序员来说这是一套开箱即用的数学计算工具能有效提升项目开发效率。1. newmat.lib 放进 VC 工程里不只是“能算”拿到这份 newmat.lib.zip 时第一反应是“又是一个被打包好的旧库”。但把压缩包翻开后才注意到里面不仅有newmat.lib这个编译产物和newmat.h主头文件还包括newmat1.cpp到newmat9.cpp的全部实现源码以及cholesky.cpp、svd.cpp、jacobi.cpp、hholder.cpp、evalue.cpp、newfft.cpp这些专门做数值算法的小模块。对用 Visual C 下手的人来说这意味着两件事一是不用为 BLAS/LAPACK 的环境挨个配置二是不必依赖别人打包那一个.lib的编译环境源码在手随时能重编出匹配自己工程 ABI 的静态库。接下来会把库的源码结构、VC 集成方式、算法模块边界和打包示例代码都过一遍每个步骤说明为什么这么设。2. 拆包后先分清源码构成再决定用哪条路径进 VS 工程2.1 三类文件别混编先建立映射解压后文件很多先按角色划分。newmat.h是主头文件newmatap.h是高级数值接口newmatio.h提供输出include.h在编译期调整 bool、STL 等兼容性设置。实现层分两组newmat1.cpp到newmat9.cpp加上bandmat.cpp、submat.cpp、myexcept.cpp是基础矩阵运算cholesky.cpp、svd.cpp、jacobi.cpp、hholder.cpp、evalue.cpp、newfft.cpp、sort.cpp、nm_misc.cpp、newmatnl.cpp是数值算法。剩下tmt*.cpp、nm_ex*.cpp、garch.cpp、nl_ex.cpp、sl_ex.cpp是测试和示例不要混进库工程。角色主要文件在 Visual Studio 工程里的放置基础类型头文件include.h, newmat.h, newmatrc.h, myexcept.hinclude 目录高级算法头文件newmatap.h, newmatnl.hinclude 目录核心实现newmat1.cpp~newmat9.cpp, bandmat.cpp, submat.cpp, myexcept.cpp编入静态库数值算法实现cholesky.cpp, svd.cpp, jacobi.cpp, hholder.cpp, evalue.cpp, newfft.cpp, nm_misc.cpp, newmatnl.cpp按需编入静态库测试与示例tmt1.cpp~tmtm.cpp, nm_ex1/2/3.cpp, garch.cpp, nl_ex.cpp, sl_ex.cpp独立可执行工程不进库这份映射是我在多个用 newmat 的老工程里验证过的。最容易踩的问题是把tmt1.cpp也加进库工程然后发现符号重复和额外的 main 入口。库工程只产出静态库测试代码单独生成 exe。2.2 Visual Studio 工程方式编出供自己用的 newmat.lib在 VS 里新建“静态库”项目目标平台选 x64。把上一节第二、三行所有实现文件拖进“源文件”头文件放“头文件”目录。配置 Debug/Release 两套都编译要特别注意语言选项设置“C/C ▸ 代码生成”里的运行时库为“多线程 DLL (/MD)”这样使用方工程无论走 Release 还是 Debug 都能链上对应的 CRT。异常模型保持/EHscnewmat 内部大量依赖异常上报越界。编译选项做一句话定位/O2开启速度优化newmat 这类数值库在 O2 下运行效率接近普通手写裸循环/MD是大多数现代应用项目的选择如果换成/MT静态 CRT 也不是不能跑但后续链接另一个也用了/MT的第三方库时会频繁报 LNK4098 警告。如果最终工程使用 VS2015 之后的版本建议把库工程和使用方工程的“平台工具集”保持在同一版本避免运行时行为不一致。2.3 命令行快速重编一次CI 里也能用没有打开 IDE 时包里的nm_i8.mak、nm_b55.mak、nm_gnu.mak、nm_ow.mak分别对应 Intel、Borland、GCC、OpenWatcom但没提供一个现成的 MSVC makefile。我一般在 CI 脚本里直接调cl和lib两条命令cl /nologo /O2 /EHsc /MD /c newmat1.cpp newmat2.cpp newmat3.cpp newmat4.cpp \ newmat5.cpp newmat6.cpp newmat7.cpp newmat8.cpp newmat9.cpp bandmat.cpp \ submat.cpp myexcept.cpp cholesky.cpp svd.cpp jacobi.cpp hholder.cpp \ evalue.cpp newfft.cpp sort.cpp nm_misc.cpp newmatnl.cpp lib /NOLOGO /OUT:newmat.lib *.obj第一行把矩阵库全部实现源码分别编译成对象文件。参数含义很直接/c代表只编译不链接免得cl因为入口函数缺失报 LNK2019/O2对应速度优化对大型矩阵乘法收益较大/MD使用动态 CRT和主工程保持一致。第二行的lib是 MSVC 自带的库管理器把上一行产出的所有.obj包进newmat.lib。整个过程不用配置 VS 工程文件之后无论是 CMake 还是 msbuild 项目都可以直接把newmat.lib当作外部依赖链接。提示如果工程同时链接别人给的一份 newmat.lib又把newmat1.cpp等源文件作为编译单元加进去链接阶段会出现复数符号定义错误。源码和 .lib 只能二选一。3. newmat 的矩阵类型与引用语义决定你该怎么写运算3.1 矩阵类型家族Matrix仍是中心但其余类型不是它的子类newmat 这里比较反直觉的一点是SymmetricMatrix、BandMat、DiagonalMatrix不是从Matrix继承的而是都从BaseMatrix这个抽象基类派生。这样做的目的是让重载运算符的返回值类型能按具体运算保持比如SymmetricMatrix SymmetricMatrix可以直接返回SymmetricMatrix如果强迫先转成Matrix再运算带宽和对称信息就浪费了。BandMat内部只存带内的lowerupper1行元素访问(i,j)时先检查是否落在带内越界时bandmat.cpp里的函数会抬一个 index error。写通用算法时参数要写成const BaseMatrix实际计算在Matrix构造时完成展开。常见的坑是直接拿const BaseMatrix调.t()、.i()看起来能用实际返回的是临时对象再用它参与运算容易把类型退化掉。所以我在团队里约定函数入参用const BaseMatrix但对入参做的最终运算必须在函数体内显式转成具体类型。运算对象类型存储内容常见误用Matrix全矩形数据行主序拿 0 当行下标SymmetricMatrix下三角存储不包括上三角写S(i,j)且i j之外的元素BandMat只存带宽内元素对带外元素赋值DiagonalMatrix单一主对角数组访问D(i,j)时i ! j3.2 运算符的临时对象机制不是所有表达式都能随意写newmat 内部做了完整的表达式临时对象复用A B C D这类长表达式能编译也能跑但如果你把同一矩阵放在赋值号两边例如A A B它能处理。真正危险的是A A A.t()这种在临时对象构造过程中复用了同一块存储的情况求转置时可能读到已经覆盖过的数据。安全示例是这样的Matrix A(2, 2); A 1.0 2.0 3.0 4.0; Matrix B A A.t(); // 安全临时对象类型明确 Matrix C; C A A.t(); // 也安全 Matrix D A * A; // 安全结果矩阵是独立分配的 label:这里B是拷贝构造C先默认构造再走赋值运算符两种情况都不共享数据。需要真正避免的是 in-place 的更新newmat 没有实现类似A A.t()这种自反语句优化对于SymmetricMatrix尤其容易出问题。安全做法是先用临时对象算出结果再赋回去。代码块的注释不用理会label:只是给断点用的一个空标签编译型代码里可以删掉。3.3 子矩阵获取columns() 返回的不是裸引用submat.cpp里实现了一套较完整的子矩阵接口A.rows(a, b)、A.columns(c, d)、A.submatrix(r1, r2, c1, c2)都返回可分派给Matrix的结果对象。这里很多人被“视图”两个字误导以为拿到的是不需要拷贝的引用。其实 newmat 的做法是返回一个持有原矩阵引用的临时对象如果把它直接赋给一个新的Matrixnewmat 会做一次完整复制。如果试图用非 const 引用去绑定这个临时对象编译器直接报错Matrix A(4, 4); Matrix M A.columns(1, 2); // 安全M 获得复制后的独立数据 // Matrix M2 A.columns(1, 2); // 错误临时对象不能绑定非 const 引用推荐做法是子矩阵赋给Matrix时就复制后续对M的任何修改不影响A。如果确实要原地更新一行用MatrixRow系列接口比如MatrixRow mr(A, 1); mr A.Row(2);这会把第二行的数据复制到第一行而不会产生新的矩阵对象。这样写语义清楚也符合 newmat 内部引用的生命周期。4. Cholesky、SVD、特征值与 FFTnewmatap 模块的调用边界4.1 Cholesky 分解构造、判断和反代求 solvecholesky.cpp对应newmatap.h里的Cholesky类。标准用法是构造时就完成分解。构造函数有两个参数第二个是可选的最小主元容忍度 epsilon一般保持默认即可。#include newmat.h #include newmatap.h SymmetricMatrix A(3); A 4.0 -2.0 -2.0 5.0 1.0 6.0; Cholesky chol(A); if (!chol.IsValid()) { cerr matrix not SPD endl; return -1; } ColumnVector b(3); b 8.0 0.0 7.0; ColumnVector x chol.i(b);在调用顺序上要提醒chol.i(b)不是先求逆矩阵再乘向量而是内部直接做 forward/back substitution得到线性方程组的解。随后取chol.cholesky()也能正常工作因为分解结果已缓存在对象内部。IsValid()的返回值基于分解过程中对角线元素与 epsilon 的比较判断所以一个理论上正定但条件数很差的矩阵可能被判无效这一点和 MATLAB 的 chol 默认行为不一样遇到边界矩阵时要有心理准备。4.2 SVD 与伪逆奇异值截断要自己处理svd.cpp提供SVD(Matrix, DiagonalMatrix, Matrix, Matrix)。它输出的是 V不是 VT。写最小二乘时我通常手动构造逆奇异值矩阵而不是调用D.i()Matrix U, V; DiagonalMatrix D; SVD(A, D, U, V); Real tol 1e-10; DiagonalMatrix Dinv(D.nrows()); for (int i 1; i D.nrows(); i) Dinv(i, i) (fabs(D(i, i)) tol) ? 1.0 / D(i, i) : 0.0; Matrix pseudoInv V * Dinv * U.t();这里Dinv(i, i)对小于 tol 的奇异值直接置 0而不是调D.i()否则 0 奇异值的倒数会变成一个异常大的数让最后结果退化。D.nrows()返回奇异值个数矩阵 U 和 V 的维度与此匹配。伪逆矩阵算出来后我再做一次pseudoInv * A的乘积检查确认结果接近单位阵避免在条件数很大的场景下自己没改干净。4.3 特征值Jacobi、Householder、EigenValues 的选用evalue.cpp、jacobi.cpp、hholder.cpp三个文件对应三套算法。实际使用中EigenValues是对称矩阵全特征值的首选接口SymmetricMatrix S(3); S 2.0 -1.0 0.0 2.0 -1.0 1.0; DiagonalMatrix D; Matrix V; EigenValues(S, D, V);D 中特征值按从大到小排序V 的每一列对应一个特征向量。如果数据量不大但需要高性能可用Jacobi做迭代如果矩阵维度较大且只求部分特征值先Householder化成三对角矩阵再对目标区间做二分和反迭代内存占用会低不少。这个选择顺序不是玄学Householder 的第一步是固定的正交变换代价是 O(n³)但之后求部分特征值只需要 O(n²) 量级的迭代整体比重启全量特征值快很多。4.4 FFT长度必须是 2 的幂否则别怪结果是 NaNnewfft.cpp里封装了FFT_Engine类。构造参数就是变换长度 n内部会按 n 构建旋转因子表。示例调用如下int n 8; FFT_Engine fft(n); ColumnVector signal(n), spectrum(n); for (int i 1; i n; i) signal(i) sin(2.0 * 3.14159265 * i / n); spectrum fft.Fft(signal);Fft返回频谱数据接口本身不复杂。复杂的是 n 的限制如果不是 2 的幂newfft.cpp内部的拆分逻辑会进入未定义分支轻则结果 NaN重则越界写内存。在调用前加一层检查把 n 向上对齐到不小于它的最小 2 的幂这是最稳妥的做法。任务类/函数限制求解对称正定线性方程组Cholesky必须对称正定最小二乘SVD需自行截断奇异值对称矩阵全部特征值EigenValues适用于对称矩阵对称矩阵三对角化Householder可用于部分特征值场景频谱分析FFT_Engine长度必须是 2 的幂5. 从 garch.cpp、nl_ex.cpp 到自己的统计模型示例代码的改造5.1 garch.cpp 是一个把矩阵求逆和循环放在一起的范例garch.cpp对多数人来说第一眼像统计模型但实际价值是演示怎样在既有项目里用 newmat 反复处理协方差矩阵和二次型。把它挪进 VC 工程照抄会先碰到两件事文件里大量用输出矩阵必须保证#include newmatio.h出现在所有输出之前另一个是文件使用Real类型在include.h里定义如果源码里直接写了double后续要全局替换。我一般把garch.cpp处理成目标函数只看它更新统计算子的部分#include newmat.h #include newmatap.h Real garch_log_likelihood(const ColumnVector y, const ColumnVector theta) { int T y.nrows(); Real omega theta(1); Real alpha theta(2); Real beta theta(3); Real h omega / (1.0 - alpha - beta); Real ll 0.0; for (int t 2; t T; t) { h omega alpha * y(t - 1) * y(t - 1) beta * h; ll -0.5 * (log(2.0 * 3.141592653589793) log(h) y(t) * y(t) / h); } return ll; }这里y是从garch.dat里读出的收益序列theta是三个待估参数。函数体本身没直接用矩阵运算但它验证了整个依赖链多变量 GARCH 的状态方程必然出现矩阵乘法和 Cholesky 求逆而这正是 newmat 的长处。zip 里同时放着nm_i10.mak、nm_m6.mak、nm_cc.mak也是提醒你这个库既可以被拆成单文件参与编译也可以交给多种 C 编译器。5.2 nl_ex.cpp 和 sl_ex.cpp不是玩具是调试工具链nl_ex.cpp演示的是NonLinearLeastSquares类把目标函数写成ReturnMatrix形式由库自己处理数值梯度和步长。这类代码看着是示例实际价值在于它提供了全局函数的调用范式当你在本地把nl_ex.cpp编译通过newmat 的异常捕获、临时矩阵管理、以及include.h里的编译器兼容开关就都已验证过。sl_ex.cpp则配合nm_misc.cpp里的Simplex函数给低维参数优化兜底。遇到目标函数梯度不稳定时先拿单纯形法跑一遍找初始值再用带梯度的优化器精修往往比直接上全局搜索收敛得更快。5.3 把示例改成可复用的封装只暴露业务接口从示例提炼规则时第一原则是不要在业务代码里散落 newmat 类型。一个最小封装是把 GARCH 的预测过程包成类newmat 类型全部留在 .cpp 内部class GarchModel { public: void setData(const std::vectordouble data); std::vectordouble predict(unsigned int horizon) const; private: ColumnVector r_; Real omega_, alpha_, beta_; }; void GarchModel::setData(const std::vectordouble data) { r_.ReSize(static_castint(data.size())); for (int i 0; i static_castint(data.size()); i) r_(i 1) data[i]; }注意r_.ReSize(n)不保证保留旧值与 STL 的resize行为不同。如果调用ReSize后想保留前几行旧数据需要先把旧值复制到临时数组。这个细节在实际更换数据源时踩过我调完ReSize后想保留前三行旧值结果全部归零。内部实现里predict再调用 5.1 的garch_log_likelihood配合优化器迭代外部业务只看到std::vector将来换掉底层库也不动接口。6. 用 tmt*.cpp 和 test_exc.cpp 验收库再谈三个旧库特有的坑6.1 先跑整套测试别急着量性能newmat 自带的tmt1.cpp到tmtm.cpp是针对不同模块的回归测试test_exc.cpp验证异常分支。把它们从库工程拆出来单独建一个 console 项目main 写在某个 tmt 文件里一次只编一个 tmt 文件逐个运行。最简单的验收命令是在 VS 开发命令行里做cl /nologo /O2 /EHsc newmat1.cpp newmat2.cpp newmat3.cpp newmat4.cpp \ newmat5.cpp newmat6.cpp newmat7.cpp newmat8.cpp newmat9.cpp bandmat.cpp \ submat.cpp myexcept.cpp cholesky.cpp svd.cpp jacobi.cpp hholder.cpp \ evalue.cpp newfft.cpp sort.cpp nm_misc.cpp newmatnl.cpp tmt1.cpp /Fe:tmt1.exe tmt1.exe正常输出末尾会出现类似 Ok 或 Tests completed 的信息。如果看到 fail多半不是库本身问题而是编译器对模板表达式的支持不同优先检查include.h顶部的条件编译宏。tmt 系列跑通后再有新的矩阵运算改动就能拿它们当基线回归集。6.2 三个旧库特有的坑绕不过去第一内存泄漏报告是性能问题不是真错误。newmat 自带 new/delete 统计控制台关闭前输出 Return from main - check for possible leaks很多初次使用者当失败处理。实际上只要连续两次运行泄漏数量不增长就可以认为正常。第二预编译的 newmat.lib 版本与主工程 CRT 不一致时异常流向会乱。压缩包里那个 newmat.lib 可能是多年前编译器编的和你当前 VS 的 C ABI 不匹配不要迷信它能直接链。拿到源码的价值就在于此重新执行第 2 节的cl和lib两条命令确保本地产物与主工程同一套工具链编译。第三USE_LAPACK宏。如果用源码编译却预定义了USE_LAPACK编译器会去找blas.h和lapack.h而包里并没有这两个头文件。检查include.h或工程预处理定义把没用的USE_LAPACK删掉很多 CI 报错都是因为这个宏从全局配置里传了下来。本文还有配套的精品资源点击获取