#matrix #condition #ndarray #norm #linalg

condest

Higham和Tisseur于2000年实现的1-范数和条件数估算器

3个不稳定版本

0.2.1 2019年1月22日
0.2.0 2019年1月22日
0.1.0 2019年1月9日

#617 in 科学

每月21次下载
用于 expm

MIT/Apache

45KB
660

Rust中的条件数和1-范数估算

此包实现了Higham和Tisseur的矩阵1-范数估算器,算法2.4,见链接文档的第7页(1190)。它允许对单个矩阵、矩阵幂和矩阵乘积进行1-范数估算,这比显式计算它们更便宜。

它使用优秀的rust-ndarray包进行矩阵存储。

示例用法

下面的示例生成一个随机矩阵a并估算其1-范数。平均来说,这给出了相当好的结果。当然,也有一些矩阵,这个算法严重低估了实际的1-范数。更多内容请参见Higham和Tisseur

condest::normest1创建一个Normest1结构体,使用它来估算1-范数,然后丢弃它。如果您想重复估算相同尺寸矩阵的1-范数,请初始化Normest1并调用normest1normest1_pownormest1_prod

重要:您需要显式链接到BLAS + LAPACK提供程序,例如openblas_src。请参阅blas-lapack-rs组织提供的说明。

extern crate openblas_src; // Need to declare `openblas_src` (or some other BLAS provider) explicitly to link to a BLAS library.

use ndarray::{
    prelude::*,
    Zip,
};
use ndarray_rand::RandomExt;
use rand::{
    SeedableRng,
};
use rand::distributions::StandardNormal;
use rand_xoshiro::Xoshiro256Plus;

fn main() {
    let t = 2;
    let n = 100;
    let itmax = 5;

    let mut rng = Xoshiro256Plus::seed_from_u64(1234);
    let distribution = StandardNormal;

    let mut a = Array::random_using((n, n), distribution, &mut rng);
    a.mapv_inplace(|x| 1.0/x);

    let estimated = condest::normest1(&a, t, itmax);
    let expected = {
        let (layout, a_slice) = if let Some(a_slice) = a.as_slice() {
            (cblas::Layout::RowMajor, a_slice)
        } else if let Some(a_slice) = a.as_slice_memory_order() {
            (cblas::Layout::ColumnMajor, a_slice)
        } else {
            panic!("Matrix not contiguous!")
        };

        unsafe {
            lapacke::dlange(
            layout,
            b'1',
            n as i32,
            n as i32,
            a_slice,
            n as i32,
        )}
    };

    assert_eq!(estimated, expected);
}

待办事项

目前,仅公开1-范数估算。实现计算条件数所需的向量,但尚未通过API访问。待解决的问题包括

  • 返回计算1-范数所需的向量;
  • 创建一个结构体,其中包含必要的临时变量,以便在不进行额外分配的情况下重复调用normest1
  • 实现额外的测试,以模拟Higham和Tisseur中的数值实验。
  • 编写一些优秀的文档。

依赖项

~5MB
~147K SLoC