IMSL_DROP_TOLERANCE, float tol (入力)
こ 階数決定 許容値 す 候補 列：
a c d
d 最初 ン よ 下 値が消去さ います そ 結果 相対条件：
2 2
d tol c .
満たさ け い そう い場合 列 a 以前に落 さ た列に線形依存す 制約 引 続 満たさ ます
フ : tol = sqrt(imsl_f_machine(3))
IMSL_SUPPLY_WORK_ARRAYS , int lwork, float work[], int liwork, int iwork[] (入力/出力) こ プ ン引数 使うこ によ ー 作業用配列 work iwork イ 位置 提供 大 問題に対し 効率が高く 断片化も避け こ
が ます maxt アク に ッ 最大数 す 次 条件が満たさ
ませ ：
lwork ≥ maxt*(m*(n+2) + n) よび liwork ≥ maxt*n
OpenMP 並列 ッ ング 使わ い場合 maxt=1
IMSL_OPTIMIZED, int *flag (出力)
最適 残差 が得 た う 示す 0-1 フ グ 1 値 最適残差 が
得 たこ 示します 最大 反復数に達した場合 0 に ます フ グ 説明
0 最大 反復数に達した
1 最適 残差 が得 た
IMSL_DUAL_SOLUTION, float **dual (出力)
対ベク w
= A
T(Ax - b)
含 長さ n 配列 こ 反復 最大数に先に達しいた 最適 い (全 成分が w ≤ 0 満たすこ い) こ もあ ます
IMSL_DUAL_SOLUTION_USER, float dual[] (出力)
ー によ 提供さ 対 た 領域 IMSL_DUAL_SOLUTION 参照 IMSL_RESIDUAL_NORM, float *rnorm (出力)
残差ベク
|| Ax - b ||
2IMSL_RETURN_USER, float x[] (出力)
ー 割 当 長さ n 配列 近似解ベク x ≥ 0 含
説明
関数 imsl_f_nonneg_least_squares 制約 x ≥ 0
|| Ax - b ||
2 最小にす こ によ Ax ≈ b 最小二乗解 求 ます 使用す ア NNLS こ Charles L. Lawson and Richard J.Hanson, Solving Least Squares Problems, SIAM Publications, Chap. 23, (1995) に説明さ います こ 関数 が採用し い チ ッ 制約 無くすこ 新しい機能 す ジ NNLS法 チ ッ に対応し いませ また 一 けが使わ 場合 も全 対成分が計算さ ます 最初 に遭遇した変数 非活性にす 通常性能が向 します 最適 解 い 方法 得 ます m n 相対 イ に制限 あ ませ
例題 1
関数 し 次 指数関数 考えます：
f(t) = c1 + c2 exp(-λ2 t) + c3 exp(-λ3 t), t ≥ 0
指数部 引数 パ ータ
2 1, 3 5
固定さ ます 係数
0, 1,2,3 cj j ンプ ータ値
i, 1,...,21
f t i非負 最小二乗法 推定します ータ し 用い 値 ti = 0.25i, i = 1, ..., 20
ここ
c1 = 1, c2 = 0.2, c3 = 0.3
#include <imsl.h>
#include <math.h>
#define M 21
#define N 3
int main() { int i;
float a[M][N], b[M], *c;
for (i = 0; i < M; i++) {
/* Generate exponential values. This model is y(t) = c_0 + c_1*exp(-t) + c_2*exp(-5*t) */
a[i][0] = 1.0;
a[i][1] = exp(-(i*0.25));
a[i][2] = exp(-(i*0.25)*5.0);
/* Compute sample values */
b[i] = a[i][0] + 0.2*a[i][1] + 0.3*a[i][2];
}
/* Solve for coefficients, constraining values to be non-negative. */
c = imsl_f_nonneg_least_squares(M, N, &a[0][0], b, 0);
/* With noise level = 0, solution should be (1, 0.2, 0.3) */
imsl_f_write_matrix("Coefficients", 1, N, c, 0);
}
出力
Coefficients
1 2 3 1.0 0.2 0.3
例題
指数関数
1 2exp(
2)
3exp(
3) ( ), 0 f t c c t c t n t t
ここ 2,3 値 例題1 同 す 関数 n(t) 標準偏差
σ
= 10-3 正規分布す ンダ イ表わします n(t) に対し ns = 10001個 ンプ ー ン 行 いました そ
結果 問題 OpenMP 並列 解 ます OpenMP 結果が正しいこ チ ックす た に
OpenMP 使わ い ープも計算します こ 残差 一致す OpenMP 並列にした
場合 逐次処理した場合 結果に差が無いこ が分 ます
#include <imsl.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>
#define M 21
#define N 3
#define NS 10001
int main() {
#define BS(i_,j_) bs[(i_)*M + (j_)]
#define X(i_,j_) x[(i_)*N + (j_)]
int thread_safe=1, seed=123457, i, *iwork, j, lwork, liwork, maxt;
float b[M], *work, sigma=1.0e-3, a[M][N], rseq[NS], rpar[NS], *bs, *x;
/* Allocate work memory for all threads that are used in the loops below. */
maxt = omp_get_max_threads();
lwork = maxt*(M*(N+2)+N);
liwork = maxt*N;
work = (float *) malloc(lwork * sizeof(float));
iwork = (int *) malloc(liwork * sizeof(int));
x = (float *) malloc(NS*N * sizeof(float));
bs = (float *) malloc(NS*M * sizeof(float));
for (i = 0; i < M; i++) { /* Generate matrix values.
This model is y(t) =
c_0 + c_1*exp(-t) + c_2*exp(-5*t) + n(t) */
a[i][0] = 1.0;
a[i][1] = exp(-(i*0.25));
a[i][2] = exp(-(i*0.25)*5.0);
}
/* Solve for coefficients, constraining values to be non-negative.
First use a sequential for loop. Then a parallel for loop.
Record the residual norms and compare them. */
imsl_random_seed_set(seed);
/* First the sequential loop.
Working memory is not included as an argument. */
for (j = 0; j < NS; j++) {
imsl_f_random_normal(M, IMSL_RETURN_USER, b, 0);
/* Add normal pdf noise at the level sigma. */
for (i=0; i<M; i++) {
b[i] = sigma*b[i] + a[i][0] + 0.2*a[i][1] + 0.3*a[i][2];
BS(j,i) = b[i];
}
imsl_f_nonneg_least_squares(M, N, &a[0][0], &BS(j,0), IMSL_RETURN_USER, &X(j,0),
IMSL_RESIDUAL_NORM, &rseq[j],
0);
}
/* Then the parallel for loop using OpenMP.
Working memory is an optional argument. This is not required but helps prevent memory fragmentation. */
/* Reset x for output for the OpenMP loop. */
for (i = 0; i < NS*N; i++) x[i] = 0.0;
#pragma omp parallel for private(j) for (j = 0; j < NS; j++) {
imsl_f_nonneg_least_squares(M, N, &a[0][0], &BS(j,0), IMSL_RETURN_USER, &X(j,0),
IMSL_RESIDUAL_NORM, &rpar[j],
IMSL_SUPPLY_WORK_ARRAYS, lwork, work, liwork, iwork, 0);
}
/* Check that residual norms agree exactly for both loops. They should because the same problems are solved - one set
sequentially and the next set in parallel. */
for (j = 0; j < NS; j++) {
/* Since the two loops solve the same set of problems, the residual norms must agree exactly. */
if (rpar[j] != rseq[j]) { thread_safe = 0;
break;
} }
if(thread_safe)
printf("imsl_f_nonneg_least_squares is thread-safe.\n");
else
printf("imsl_f_nonneg_least_squares is not thread-safe.\n");
system("pause");
}
出力
imsl_f_nonneg_least_squares is thread-safe.
警告 ー
IMSL_MAX_NNLS_ITER_REACHED 反復数が最大値に達した 最良 解が返さ ます
“itmax” = # が使わ ました よ 大 値にす
計算が完了す もし ませ