与给定向量形成正交基的矩阵答案

【问题标题】：matrix that forms an orthogonal basis with a given vector与给定向量形成正交基的矩阵
【发布时间】：2011-03-09 13:52:35
【问题描述】：

一道线性代数题；

给定一个 k 变量范数向量 u（即 u : ||u||_2=1）你如何构造 \Gamma_u，任意 k*(k-1) 单位向量矩阵使得 (u,\Gamma_u) 形成一个正交基？

我的意思是：从计算的角度来看：你用什么算法来构造这样的矩阵？

提前致谢，

【问题讨论】：

标签： algorithm math geometry linear-algebra

【解决方案1】：

天真的方法是应用 u_0 和 k-1 个随机生成的向量的 Gram Schmidt 正交化。如果在某个时候 GS 算法生成一个零向量，那么你有一个线性依赖，在这种情况下再次随机选择向量。

但是，这种方法不稳定，向量表示中的小数值误差会被放大。然而，这个算法存在一个稳定的修改：

让a_1 = u, a_2,...a_k 成为随机选择的向量

for i = 1 to k do 
        vi = ai
end for 

for i = 1 to k do
    rii = |vi| 
    qi = vi/rii
    for j = i + 1 to k do
       rij =<qi,vj>
       vj =vj −rij*qi 
    end for
end for

结果向量v1,...vk 将是矩阵的列，v1 = u。如果在某个时刻vj 变为零，请选择一个新向量aj 并重新开始。请注意，如果向量 a2,..,ak 是随机选择的，这种情况发生的概率可以忽略不计。

【讨论】：

谢谢你：我回家试试，让你知道，Best。
您可以不从随机基数开始。 1）检查起始向量是否是规范基数中一的倍数。 2）如果是，则选择规范基中的其他n-1个向量。 3）如果不是，检查犬基的任何 n-1 个向量。 /// 没有选择错误基地的可能性。
@belisarius，你仍然可能会得到线性依赖。考虑简单的情况 u = (1,1,0)，并且 e1=(1,0,0)，e2=(0,1,0)，e3=(0,0,1)，那么 u 不是任何一个的线性倍数，但选择 e1 和 e2 作为其他两个向量将导致线性组合。一般来说，检测线性组合并不便宜，因此使用新的随机向量重新启动 Gram-schmidt 方法可能更便宜，因为这是一个非常罕见的事件
@user189035：我不这么认为，到底是什么问题？

【解决方案2】：

您可以使用 Householder 矩阵来执行此操作。查看例如http://en.wikipedia.org/wiki/Householder_reflection 和http://en.wikipedia.org/wiki/QR_decomposition

可以找到一个 Householder 矩阵 Q 以便 Q*u = e_1 （其中e_k 是除第 k 个位置的 1 外全为 0 的向量）那么如果f_k = Q*e_k，则f_k和f_1 = u形成一个正交基。（由于Q*Q = I，Q是正交的。）

所有这些关于矩阵的讨论可能会让人觉得例行公事会很贵，但事实并非如此。例如这个 C 函数，给定一个长度为 1 的向量，返回一个具有所需基的数组按列顺序，即第 i 个向量的第 j 个分量保存在 b[j+dim*i]

   double*  make_basis( int dim, const double* v)
{
    double* B = calloc( dim*dim, sizeof * B);
    double* h = calloc( dim, sizeof *h);
    double  f, s, d;
    int i, j;

    /* compute Householder vector and factor */
    memcpy( h, v, dim*sizeof *h);
    s = ( v[0] > 0.0) ? 1.0 : -1.0;
    h[0] += s;  
    f = s/(s+v[0]);

    /* compute basis */
    memcpy( B, v, dim * sizeof *v); /* first one is v */
    /* others by applying Householder matrix */
    for( i=1; i<dim; ++i)
    {   d = f*h[i];
        for( j=0; j<dim; ++j)
        {   B[dim*i+j] = (i==j) - d*h[j];
        }
    }
    free( h);
    return B;
}

【讨论】：