6

给定一堆任意向量（存储在矩阵 A 中）和半径 r，我想找到落在半径为 r 的球体内的那些向量的所有整数值线性组合。然后我会将必要的坐标存储在矩阵 V 中。因此，例如，如果线性组合

K=[0; 1; 0]

降落在我的球内，例如

if norm(A*K) <= r then
   V(:,1)=K
end

等等

A 中的向量肯定是给定格的最简单的基础，最大向量的长度为 1。不确定这是否以任何有用的方式限制向量，但我怀疑它可能。- 他们不会像不太理想的基础那样有相似的方向。

我已经尝试了一些方法，但没有一个看起来特别令人满意。我似乎找不到一个很好的模式来遍历格子。

我目前的方法是从中间开始（即所有 0 的线性组合）并一一通过必要的坐标。它涉及存储一堆额外的向量来跟踪，所以我可以遍历坐标的所有八分圆（在 3D 情况下）并一个接一个地找到它们。这个实现看起来非常复杂并且不是很灵活（特别是它似乎不容易推广到任意数量的维度 - 尽管这对于当前目的并不是绝对必要的，但它会是一个不错的选择）

有没有一种很好的*方法来找到所有必需的点？

（*理想情况下既高效又优雅**。如果真的有必要，在球体之外有几个额外的点并不重要，但最好不要那么多。我确实需要球体内的所有向量。-如果它有很大的不同，我对 3D 案例最感兴趣。

**我很确定我目前的实现都不是。）

我发现了类似的问题：

在任意坐标周围找到半径为 r 的球体中的所有点- 这实际上是比我正在寻找的更普遍的情况。我只处理周期性格子，我的球体始终以 0 为中心，与格子上的一个点重合。但是我没有一个点列表，而是一个向量矩阵，我可以用它来生成所有点。

如何有效地枚举 n 维网格中的所有球体点- 完全规则的超立方晶格和曼哈顿距离的情况。我正在寻找完全任意的格子和欧几里得距离（或者，为了提高效率，显然是平方）。

4

4 回答 4

1

我找到了一种让我现在更快乐的方法。可能还有改进的余地，所以如果你有更好的方法，或者发现这段代码有错误，请一定分享。虽然这是我现在所拥有的：（全部用 SciLab 编写）

第 1 步：找出由与晶格向量的轴对齐的有界 n 平行四边形所定义的最大范围。感谢 ElKamina 的含糊建议以及对我在 math.se 上的另一个问题的回复：https ://math.stackexchange.com/a/1230160/49989

    function I=findMaxComponents(A,r) //given a matrix A of lattice basis vectors
                                      //and a sphere radius r,
                                      //find the corners of the bounding parallelotope
                                      //built from the lattice, and store it in I.
        [dims,vecs]=size(A); //figure out how many vectors there are in A (and, unnecessarily, how long they are)
        U=eye(vecs,vecs); //builds matching unit matrix
        iATA=pinv(A'*A); //finds the (pseudo-)inverse of A^T A
        iAT=pinv(A'); //finds the (pseudo-)inverse of A^T
        I=[]; //initializes I as an empty vector
        for i=1:vecs                           //for each lattice vector,
            t=r*(iATA*U(:,i))/norm(iAT*U(:,i)) //find the maximum component such that
                                               //it fits in the bounding n-parallelotope
                                               //of a (n-1)-sphere of radius r
            I=[I,t(i)];                        //and append it to I
        end
        I=[-I;I]; //also append the minima (by symmetry, the negative maxima)
    endfunction

在我的问题中，我只要求提供一般基础，即对于 n 维，一组 n 任意但线性独立的向量。上面的代码，凭借使用伪逆，适用于任意形状的矩阵，同样，Scilab 的 " A'" 返回共轭转置，而不仅仅是转置，A因此它同样适用于复杂矩阵。

在最后一步中，我放置了相应的最小组件。

以 A 为例，这在 Scilab 的控制台中为我提供了以下信息：

     A  =
     
        0.9701425  - 0.2425356    0.
        0.2425356    0.4850713    0.7276069
        0.2425356    0.7276069  - 0.2425356
    
    r=3;

    I=findMaxComponents(A,r)
    
     I  =
     
      - 2.9494438  - 3.4186986  - 4.0826424
        2.9494438    3.4186986    4.0826424
    
     I=int(I)
     
     I  =
     
      - 2.  - 3.  - 4.
        2.    3.    4.

找到的值findMaxComponents是每个晶格向量的最大可能系数，使得存在与该系数的线性组合仍然落在球体上。由于我正在寻找具有整数系数的最大此类组合，因此我可以安全地删除小数点后的部分以获得最大合理的整数范围。因此，对于给定的矩阵A，我必须在第一个组件中从 -2 到 2，在第二个组件中从 -3 到 3，在第三个组件中从 -4 到 4，我肯定会落在所有点上在球体内部（加上多余的额外点，但重要的是内部的每个有效点）接下来：

第二步：利用以上信息，生成所有候选组合。

    function K=findAllCombinations(I) //takes a matrix of the form produced by
                                      //findMaxComponents() and returns a matrix
                                      //which lists all the integer linear combinations
                                      //in the respective ranges.
        v=I(1,:); //starting from the minimal vector
        K=[];
        next=1; //keeps track of what component to advance next
        changed=%F; //keeps track of whether to add the vector to the output
        
        while or(v~=I(2,:)) //as long as not all components of v match all components of the maximum vector
            if v <= I(2,:) then //if each current component is smaller than each largest possible component
                if ~changed then
                    K=[K;v]; //store the vector and
                end
                v(next)=v(next)+1; //advance the component by 1
                next=1; //also reset next to 1
                changed=%F;
            else
                v(1:next)=I(1,1:next); //reset all components smaller than or equal to the current one and
                next=next+1; //advance the next larger component next time
                changed=%T;
            end
        end
        K=[K;I(2,:)]'; //while loop ends a single iteration early so add the maximal vector too
                       //also transpose K to fit better with the other functions
    endfunction

所以现在我有了这个，剩下的就是检查给定的组合是否确实位于球体内部或外部。我要做的就是：

步骤 3：过滤组合以找到实际有效的格点

    function points=generatePoints(A,K,r)
        possiblePoints=A*K; //explicitly generates all the possible points
        points=[];
        for i=possiblePoints
            if i'*i<=r*r then //filter those that are too far from the origin
                points=[points i];
            end
        end
    endfunction

我得到了所有真正适合半径 r 的球体的组合。

对于上面的示例，输出相当长：在半径为 3 的球体的最初 315 个可能点中，我得到了剩余的 163 个点。

前4个是：（每列一个）

  - 0.2425356    0.2425356    1.2126781  - 0.9701425
  - 2.4253563  - 2.6678919  - 2.4253563  - 2.4253563
    1.6977494    0.           0.2425356    0.4850713

所以剩下的工作就是优化。据推测，其中一些循环可以更快，特别是随着维度数量的增加，我必须生成大量必须丢弃的点，所以也许有比采用边界-的n-parallelotope更好的方法n-1球体作为起点。

于 2015-04-15T22:54:48.000 回答
1

顺便说一句，在没有证明任何断言的情况下，我认为 1）如果向量集不是最大等级，那么解决方案的数量是无限的；2）如果该集合是最大秩的，则向量生成的线性变换的图像是目标空间的一个子空间（例如，平面），它在一个低维球体中与球体相交；3）因此可以将问题简化为1-1线性变换（k维空间上的kxk矩阵）；4）由于矩阵是可逆的，您可以将球体“拉回”到包含晶格点的空间中的椭球体，作为奖励，您可以获得椭球体的良好几何描述（主轴定理）；5）您的问题现在正好成为确定椭球内晶格点的问题之一。

后一个问题与 Gauss 考虑的一个老问题（计算椭圆内的晶格点）有关，他得出了一个很好的近似值。确定椭圆（oid）内的晶格点可能不是一个整洁的问题，但它可能可以一次减少一个维度（椭圆体和平面的横截面是另一个椭圆体）。

于 2015-04-16T00:42:18.467 回答
0

让我们将 K 表示为 X。

问题可以表示为：

(a11x1 + a12x2..)^2 + (a21x1 + a22x2..)^2 ... < r^2

(x1,x2,...) 不会形成球体。

于 2015-04-09T20:00:24.840 回答
0

这可以通过维度上的递归来完成——选择一个晶格超平面方向并索引所有与 r 半径球相交的超平面。每个这样的超平面的球交点本身就是一个球，在一个较低的维度上。重复。这是 Octave 中的调用函数代码：

function lat_points(lat_bas_mx,rr)

% **globals for hyperplane lattice point recursive function**
clear global; % this seems necessary/important between runs of this function
global MLB;
global NN_hat;
global NN_len;
global INP; % matrix of interior points, each point(vector) a column vector
global ctr; % integer counter, for keeping track of lattice point vectors added 
    %   in the pre-allocated INP matrix; will finish iteration with actual # of points found

ctr = 0;    % counts number of ball-interior lattice points found
MLB = lat_bas_mx;
ndim = size(MLB)(1);

% **create hyperplane normal vectors for recursion step**
% given full-rank lattice basis matrix MLB (each vector in lattice basis a column),
%   form set of normal vectors between successive, nested lattice hyperplanes;
%   store them as columnar unit normal vectors in NN_hat matrix and their lengths in NN_len vector
NN_hat = [];
for jj=1:ndim-1
    tmp_mx = MLB(:,jj+1:ndim);
    tmp_mx = [NN_hat(:,1:jj-1),tmp_mx];
    NN_hat(:,jj) = null(tmp_mx'); % null space of transpose = orthogonal to columns
    tmp_len = norm(NN_hat(:,jj));
    NN_hat(:,jj) = NN_hat(:,jj)/tmp_len;
    NN_len(jj) = dot(MLB(:,jj),NN_hat(:,jj));
    if (NN_len(jj)<0) % NN_hat(:,jj) and MLB(:,jj) must have positive dot product 
        %   for cutting hyperplane indexing to work correctly
        NN_hat(:,jj) = -NN_hat(:,jj);
        NN_len(jj) = -NN_len(jj);
    endif
endfor
NN_len(ndim) = norm(MLB(:,ndim));
NN_hat(:,ndim) = MLB(:,ndim)/NN_len(ndim); % the lowest recursion level normal
    %   is just the last lattice basis vector

% **estimate number of interior lattice points, and pre-allocate memory for INP**
vol_ppl = prod(NN_len); % the volume of the ndim dimensional lattice paralellepiped 
    %   is just the product of the NN_len's (they amount to the nested altitudes 
    %   of hyperplane "paralellepipeds")
vol_bll = exp( (ndim/2)*log(pi) + ndim*log(rr) - gammaln(ndim/2+1) ); % volume of ndim ball, radius rr
est_num_pts = ceil(vol_bll/vol_ppl); % estimated number of lattice points in the ball
err_fac = 1.1; % error factor for memory pre-allocation--assume max of err_fac*est_num_pts columns required in INP
INP = zeros(ndim,ceil(err_fac*est_num_pts));

% **call the (recursive) function**
% for output, global variable INP (matrix of interior points)
%   stores each valid lattice point (as a column vector)
clp = zeros(ndim,1); % confirmed lattice point (start at origin)
bpt = zeros(ndim,1); % point at center of ball (initially, at origin)
rd = 1; % initial recursion depth must always be 1
hyp_fun(clp,bpt,rr,ndim,rd);

printf("%i lattice points found\n",ctr);
INP = INP(:,1:ctr); % trim excess zeros from pre-allocation (if any)

endfunction

关于 NN_len(jj)*NN_hat(:,jj) 向量——它们可以被视为由格基 MLB 中的向量形成的 ndim 维“平行六面体”中的连续（嵌套）高度。晶格基平行六面体的体积只是 prod(NN_len)——为了快速估计内部晶格点的数量，将半径为 rr 的 ndim 球的体积除以 prod(NN_len)。这是递归函数代码：

function hyp_fun(clp,bpt,rr,ndim,rd)

%{
    clp = the lattice point we're entering this lattice hyperplane with
    bpt = location of center of ball in this hyperplane
    rr = radius of ball
    rd = recrusion depth--from 1 to ndim
%}

    global MLB;
    global NN_hat;
    global NN_len;
    global INP;
    global ctr;

    % hyperplane intersection detection step
    nml_hat = NN_hat(:,rd);
    nh_comp = dot(clp-bpt,nml_hat); 
    ix_hi = floor((rr-nh_comp)/NN_len(rd));
    ix_lo = ceil((-rr-nh_comp)/NN_len(rd));
    if (ix_hi<ix_lo)
        return % no hyperplane intersections detected w/ ball; 
                %   get out of this recursion level
    endif
    hp_ix = [ix_lo:ix_hi]; % indices are created wrt the received reference point
    hp_ln = length(hp_ix);

    % loop through detected hyperplanes (updated)
    if (rd<ndim)

        bpt_new_mx = bpt*ones(1,hp_ln) + NN_len(rd)*nml_hat*hp_ix; % an ndim by length(hp_ix) matrix
        clp_new_mx = clp*ones(1,hp_ln) + MLB(:,rd)*hp_ix; % an ndim by length(hp_ix) matrix
        dd_vec = nh_comp + NN_len(rd)*hp_ix; % a length(hp_ix) row vector
        rr_new_vec = sqrt(rr^2-dd_vec.^2);

        for jj=1:hp_ln
            hyp_fun(clp_new_mx(:,jj),bpt_new_mx(:,jj),rr_new_vec(jj),ndim,rd+1);
        endfor

    else  % rd=ndim--so at deepest level of recursion; record the points on the given 1-dim 
          %   "lattice line" that are inside the ball

        INP(:,ctr+1:ctr+hp_ln) = clp + MLB(:,rd)*hp_ix;
        ctr += hp_ln;
        return

    endif

endfunction

这里面有一些 Octave-y/Matlab-y 的东西，但大多数应该很容易理解；M(:,jj) 引用矩阵 M 的 jj 列；tic ' 表示转置；[AB] 连接矩阵 A 和 B；A=[] 声明一个空矩阵。

从原始答案更新/更好地优化：

  • “向量化”递归函数中的代码，以避免大多数“for”循环（那些将其减慢了约 10 倍；但现在的代码有点难以理解）
  • 为 INP 矩阵内部点预分配内存（这将其加速了另一个数量级；在此之前，Octave 必须为每次调用最内层递归级别调整 INP 矩阵的大小——对于大型矩阵/数组这真的可以减慢速度）

因为这个例程是一个项目的一部分，所以我也用 Python 编写了它。从非正式测试来看，Python 版本比这个（Octave）版本快 2-3 倍。

作为参考，这里是这个答案的原始帖子中旧的、慢得多的代码：

    % (OLD slower code, using for loops, and constantly resizing
    %   the INP matrix) loop through detected hyperplanes
    if (rd<ndim)
        for jj=1:length(hp_ix)
            bpt_new = bpt + hp_ix(jj)*NN_len(rd)*nml_hat;
            clp_new = clp + hp_ix(jj)*MLB(:,rd);
            dd = nh_comp + hp_ix(jj)*NN_len(rd);
            rr_new = sqrt(rr^2-dd^2);
            hyp_fun(clp_new,bpt_new,rr_new,ndim,rd+1);
        endfor
    else  % rd=ndim--so at deepest level of recursion; record the points on the given 1-dim
          %   "lattice line" that are inside the ball
        for jj=1:length(hp_ix)
            clp_new = clp + hp_ix(jj)*MLB(:,rd);
            INP = [INP clp_new];
        endfor
        return
    endif
于 2021-10-31T18:26:10.537 回答