Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]


Groups > comp.programming > #16283

Re: Another little puzzle

From "Dmitry A. Kazakov" <mailbox@dmitry-kazakov.de>
Newsgroups comp.programming
Subject Re: Another little puzzle
Date 2023-01-10 09:06 +0100
Organization Aioe.org NNTP Server
Message-ID <tpj6ac$56n$1@gioia.aioe.org> (permalink)
References (11 earlier) <877cxwwyhc.fsf@bsb.me.uk> <tpfcn6$1mp5$1@gioia.aioe.org> <87v8lgv17t.fsf@bsb.me.uk> <tpgpu1$1asl$1@gioia.aioe.org> <87k01vv40e.fsf@bsb.me.uk>

Show all headers | View raw


On 2023-01-09 21:37, Ben Bacarisse wrote:
> "Dmitry A. Kazakov" <mailbox@dmitry-kazakov.de> writes:
> 
>> On 2023-01-09 04:25, Ben Bacarisse wrote:
>>> "Dmitry A. Kazakov" <mailbox@dmitry-kazakov.de> writes:
>>>
>>>> On 2023-01-08 21:41, Ben Bacarisse wrote:
>>>>
>>>>> But
>>>>> the problem is /finding/ a specific average -- the point (or angle) that
>>>>> minimises the sum of squares of the distances (or angles) from that
>>>>> average point (or angle).
>>>>
>>>> That is not a programming problem. It is a problem space's one. You
>>>> must go to the corresponding part of science or engineering or
>>>> economics etc and formulate it there in terms of that space.
>>> I don't follow.  It's a mathematical problem for which an algorithm is
>>> sought.
>>
>> Firstly, this is comp.programming group.
> 
> Do you consider discussion algorithms off topic here?  Seriously?

No, I consider off topic discussing mathematical problems /= algorithms. 
You did not discuss any algorithms so far.

>> Secondly, if algorithm is sought then there must be a mathematical
>> formulation of what a numeric solution is sought for.
> 
> There needs to be a problem statement that can be discussed and refined.
> The original post about circular averages has led to two different
> refinements -- the vector average (projected onto the unit circle) and
> the arc-length average.  Both generalise to more dimensions and I
> happened to mention in passing that the great-circle average would be a
> natural fit in some situations.

I missed generalization to 3D. Sorry.
[You kept on writing about unit *circle* and an angle, where it actually 
meant unit sphere and polar angles/spherical coordinates]

> Tim mentioned (again in passing) that it looks hard.  You seemed to
> suggest it was not, but I am not sure you knew what the problem was at
> that time.
> 
> All of this seems to be to be exactly what comp.programming should be
> about.

Differential geometry and spherical trigonometry are certainly off topic.

>> (Finding a minimum is an optimization problem, there are thousands of
>> algorithms already)
> 
> Programmers should know that specifications are not algorithm designs.

Who said they are?

>>> I don't see any connection to some science or engineering
>>> space.  If it /came/ from some scientific study, it has already been
>>> turned into what the programmer needs come up with.
>>>
>>>> Average is just such formulation of finding a representative member
>>>> from some set having some specific properties. E.g. least squares in
>>>> physics usually has the meaning of least energy etc. So, it goes in
>>>> this order (not in reverse):
>>>>
>>>> 1. You need the least energy representative.
>>>> 2. An average gives you that.
>>>> 3. You measure the inputs.
>>>>
>>>> After that you ask yourself given these measures how to compute the
>>>> average? And only now it becomes a programming issue.
>>> So if I asked you for code to calculate the arithmetic mean of an array
>>> of numbers you would ask for all this first, declaring it not yet a
>>> programming issue?
>>
>> Then I would ask you, if this is a solution, what was the problem?
> 
> This is a solution to the simpler problem.

And what was that the problem? [Adjective adds/clarifies nothing]

> But we are not in that situation.  A few knowledgeable and experienced
> people are discussing a programming problem.

I saw no discussion on programming so either.

> You are welcome to join
> in, but the best way for you to do that is to ask whatever questions you
> think would help you understand what's being discussed.  So far, your
> replies have alternated between "it's trivial" and "it's unstated".

Yes, because you wrote about minimum of arc length on a *circle*. And 
that is trivial.

>>>>> You say you need a formula, so I'll try...  Let P_n be a collection of n
>>>>> unit vectors specifying n points on a unit sphere.  Find the unit vector
>>>>> A that minimises
>>>>>      Sum_{i=1,n} ( arctan( |A x P_n|  /  A . P_n ) )^2
>>>>> (arctan is the "all quadrant" version that is often called atan2 in
>>>>> programming languages.)
>>>>
>>>> 1. atan2 has two arguments (the adjacent and the opposite legs);
>>> Obviously I can't guess how much I can safely assume so I expected there
>>> might be questions.  arctan(p / q) is usually coded as atan2(p, q).
>>>
>>>> 2. What is "A";
>>> The unit vector that is being sought -- the result of the algorithm.  It
>>> represents a point on the unit sphere that is the desired average of the
>>> input points.
>>>
>>>> 3. What operations(?) "x" and "." denote?
>>> x denotes the cross product, and . the dot product.
>>> Do you agree that this is not a trivial problem?  If not, please give
>>> the trivial algorithm!
>>
>> If I correctly understood your explanation (and with a few obvious
>> corrections):
>>
>>     Sum atan2 (|A x Pi|, A · Pi)  --> min
>>    i=1..n                       {A | A on the circle}
> 
> No.  A is a vector denoting a point on the unit sphere.  And you have
> missed out the power of two.  Were these changes what you call "obvious
> corrections"?
> 
>> |A x Pi| = |A| |Pi| |sin θi|  (θi is the angle between the vectors A and Pi)
>>
>> A · Pi = |A| |Pi| cos θi
>>
>> atan2 (|A x Pi|, A · Pi) = atan2 (|sin θi|, cos θi) = θi'
>>
>> Where θi' is θi with lower semicircle mapped to the upper one (for
>> whatever reason).
>>
>> (If mapping was unintended you must have taken the direction of the
>> cross product for the sign for |A x Pi|)
>>
>> On a circle varying the angle between A and (R, 0), where R is the
>> radius, you find that the mean of the angles between Pi and (R, 0)
>> yields the minimum.
> 
> This does not even solve the related case of finding the average point
> on the unit circle.

It elementary does for both formula you gave (after corrections) and the 
circumference:

    Sum Length_Of_Arc (Pi, A)^2 =
   i=1..n

=  Sum (R Angle (Pi, A))^2 =
   i=1..n

=  R^2 Sum ((Angle (Pi, (R, 0)) - Angle (A, (R, 0)))^2
       i=1..n

(R, 0) denotes vector x = R, y = 0. Let θi denote Angle (Pi, (R, 0)) and 
α do Angle (A, (R, 0)) then

=  R^2 Sum (θi - α)^2 --> min
       i=1..n

Derivative d/dα:

   2 R^2 Sum (θi - α) = 0 <=> (R /= 0)
        i=1..n

<=> Sum θi - nα = 0 <=>
    i=1..n

<=> α = Sum θi / n = mean ({θi})
        i=1..n

The angle between A and (R, 0) is the mean of angles of points.

I can't help you with spherical trigonometry, which is an utter off 
topic anyway. I do not know if there is a direct solution of least 
squares. For two points it is like for a circle. For three points it 
would be the circumcenter of the spherical triangle.

However even if the direct solution existed, it could be numerically 
unstable as many formulae in spherical geometry are. I used to process 
remote sensing geological data which involved miles long chains of 
cosines and sines. That stuff was no fun.

Anyway, there are lengthy papers on spherical averages introducing 
iterative solutions, some with source code for the algorithms. That is 
the right way to compute it, if you really wanted...

[Again, it is not programming. In case of spherical geometry a 
programmer will find an expert mathematician and who will have the 
problem properly stated. The next stop is by an expert applied 
mathematician who will outline a workable numeric solution. Then the 
programmer can start.]

-- 
Regards,
Dmitry A. Kazakov
http://www.dmitry-kazakov.de

Back to comp.programming | Previous | NextPrevious in thread | Next in thread | Find similar


Thread

Re: Another little puzzle Ben Bacarisse <ben.usenet@bsb.me.uk> - 2022-12-30 01:00 +0000
  Re: Another little puzzle Tim Rentsch <tr.17687@z991.linuxsc.com> - 2022-12-29 23:25 -0800
    Re: Another little puzzle Ben Bacarisse <ben.usenet@bsb.me.uk> - 2022-12-30 14:04 +0000
      Re: Another little puzzle Ben Bacarisse <ben.usenet@bsb.me.uk> - 2022-12-31 00:24 +0000
        Re: Another little puzzle Tim Rentsch <tr.17687@z991.linuxsc.com> - 2022-12-31 06:42 -0800
          Re: Another little puzzle "Dmitry A. Kazakov" <mailbox@dmitry-kazakov.de> - 2022-12-31 17:04 +0100
            Re: Another little puzzle Ben Bacarisse <ben.usenet@bsb.me.uk> - 2023-01-01 01:24 +0000
            Re: Another little puzzle Tim Rentsch <tr.17687@z991.linuxsc.com> - 2023-01-08 07:45 -0800
              Re: Another little puzzle "Dmitry A. Kazakov" <mailbox@dmitry-kazakov.de> - 2023-01-08 17:17 +0100
                Re: Another little puzzle Ben Bacarisse <ben.usenet@bsb.me.uk> - 2023-01-08 20:41 +0000
                Re: Another little puzzle Richard Heathfield <rjh@cpax.org.uk> - 2023-01-08 21:14 +0000
                Re: Another little puzzle "Dmitry A. Kazakov" <mailbox@dmitry-kazakov.de> - 2023-01-08 22:31 +0100
                Re: Another little puzzle Ben Bacarisse <ben.usenet@bsb.me.uk> - 2023-01-09 03:25 +0000
                Re: Another little puzzle "Dmitry A. Kazakov" <mailbox@dmitry-kazakov.de> - 2023-01-09 11:22 +0100
                Re: Another little puzzle Ben Bacarisse <ben.usenet@bsb.me.uk> - 2023-01-09 20:37 +0000
                Re: Another little puzzle "Dmitry A. Kazakov" <mailbox@dmitry-kazakov.de> - 2023-01-10 09:06 +0100
                Re: Another little puzzle Ben Bacarisse <ben.usenet@bsb.me.uk> - 2023-01-11 02:41 +0000
                Re: Another little puzzle "Dmitry A. Kazakov" <mailbox@dmitry-kazakov.de> - 2023-01-11 10:01 +0100
                Re: Another little puzzle Ben Bacarisse <ben.usenet@bsb.me.uk> - 2023-01-12 01:00 +0000
                Re: Another little puzzle Tim Rentsch <tr.17687@z991.linuxsc.com> - 2023-01-10 23:36 -0800
                Re: Another little puzzle Y A <ya00000100000@yahoo.com> - 2023-01-11 02:39 -0800
          Re: Another little puzzle Ben Bacarisse <ben.usenet@bsb.me.uk> - 2023-01-01 01:10 +0000
            Re: Another little puzzle Tim Rentsch <tr.17687@z991.linuxsc.com> - 2023-01-08 07:17 -0800
              Re: Another little puzzle Ben Bacarisse <ben.usenet@bsb.me.uk> - 2023-01-08 19:43 +0000
                Re: Another little puzzle Ben Bacarisse <ben.usenet@bsb.me.uk> - 2023-01-08 19:59 +0000
                Re: Another little puzzle Tim Rentsch <tr.17687@z991.linuxsc.com> - 2023-01-10 23:25 -0800
      Re: Another little puzzle Tim Rentsch <tr.17687@z991.linuxsc.com> - 2022-12-31 06:20 -0800
        Re: Another little puzzle Augǝl <angel0000000001000000000000@mail.ee> - 2022-12-31 10:23 -0800

csiph-web