Excepting the question of accuracy (about which more below),
premultiplication would seem to be an apt candidate for inclusion on the
to-do list because
0. The speedup is considerable;
1. The space improvement is even more considerable, because for an mxn
matrix only a much smaller nxn product need be computed, assuming the
transpose is done implicitly during the premultiplication;
2. Using %. is hazardous anyway if y is near to having dependent
columns, so that numeric errors then might not be so important.
Just to see what the difference in error was I did a small test:
x =. (1e6 3 ?...@$ 0) mp 3 3 $ 1 0 1 0 1 1 0 0 0.000000001
y =. 1e6 1 ?...@$ 0
0j15 ": +/ (* +) x - y mp x %. y
146020.842352243870000 145878.975589646230000 417135.702066504280000
0j15 ": +/ (* +) x - y mp (%.(|:y) mp y) mp (|:y) mp x
146020.842352243870000 145878.975589646030000 417135.702066503870000
the results of the two methods are nearly identical and seem just as
often to give a smaller error to the premultiplication method as to the
QR method. I don't claim that this test proves anything other than that
the difference is small in some cases.
Does anyone here know what the worst-case operands are for the
premultiplication method? or for the QR method?
Henry Rich
On 10/10/2010 10:57 AM, Roger Hui wrote:
> When there is general code that works it takes special
> code to exploit specific sets of arguments, and special
> effort to recognize that special code would be useful.
> It appears that no special effort was taken in this case.
>
> For x%.y, the general algorithm does not require
> premultiplication by the transpose of y. Many textbooks
> say that for a tall matrix y the least-squares solution is
> (%.(|:y) mp y) mp (|:y) mp x ("premultiplication by the transpose
> of y"), but there are other, numerically more preferable
> and/or more efficient, ways. J uses the QR decomposition
> http://www.jsoftware.com/jwiki/Essays/QR%20Decomposition
> because algorithmically it is more interesting.
>
>
>
> ----- Original Message -----
> From: Henry Rich<[email protected]>
> Date: Saturday, October 9, 2010 19:34
> Subject: Re: [Jprogramming] Puzzle: line-circle/sphere intersection
> To: Programming forum<[email protected]>
>
>> I don't think I agree with your timings; what happens at rank 0
>> seems
>> unimportant.
>>
>> However, timings at full rank give the same conclusions:
>>
>> l =. 1e6 3 ?...@$ 0
>> r =. 1e6 1 ?...@$ 0
>> r %. l
>> 0.300885
>> 0.300854
>> 0.298569
>> ts 'r %. l'
>> 0.617212 1.25832e8
>> mp =: +/ . *
>> ts '(t mp r) %. (t =. |:l) mp l'
>> 0.127838 3.356e7
>> (t mp r) %. (t =. |:l) mp l
>> 0.300885
>> 0.300854
>> 0.298569
>>
>> So Roger: if you're not premultiplying by the transpose, why
>> not? (This
>> is not my area of expertise so I won't be surprised if there is
>> a good
>> reason)
>>
>> Henry Rich
>>
>> On 10/9/2010 9:02 PM, Marshall Lochbaum wrote:
>>> Bit of a problem with the matrix division thing...
>>>
>>> l=.10*<:+:1e6 ?...@$ 0
>>> r=.10*<:+:1e6 ?...@$ 0
>>> 6!:2 'l %."0 r'
>>> 3.09424
>>> 6!:2 'l (+/ .*~ %.)"0 r'
>>> 0.481168
>>> 6!:2 'l %"0 r'
>>> 0.1251
>>>
>>> It appears that multiplying by the inverse is much faster than
>> doing matrix
>>> division, which certainly shouldn't be true.
>>> It probably also explains the speed disparity between %. and
>> +/@:* -...@% ]
>>> +/@:* ] because the latter uses % in place of %. .
>>>
>>> Marshall
>>>
>>> -----Original Message-----
>>> From: [email protected]
>>> [mailto:[email protected]] On Behalf Of Henry Rich
>>> Sent: Saturday, October 09, 2010 5:22 PM
>>> To: Programming forum
>>> Subject: Re: [Jprogramming] Puzzle: line-circle/sphere intersection
>>>
>>> It seems odd to me that
>>>
>>> +/@:* -...@% ] +/@:* ]
>>>
>>> would be faster than
>>>
>>> %.
>>>
>>> for lists. I thought the first step of %. for non-
>> square arrays would
>>> be multiplying by the transpose of y (just as in the longer
>> version),> followed by taking the inverse (now of a 1x1 matrix,
>> which should be quick).
>>>
>>> Henry Rich
>>>
>>> On 10/9/2010 12:40 PM, R.E. Boss wrote:
>>>> Nice solution.
>>>> I had forgotten the projecting properties of %.
>>>>
>>>> Performance improvement is possible by substituting %. by
>> +/@:* -...@% ]
>>>> +/@:* ]
>>>>
>>>> 'C S E'=:10000 ,@#"(0 1) 0 0 0,
>> 1 2 3 ,: _2 8 1
>>>> R=:10000
>>>>
>>>> rank'C (R>:+/&.:*:)@:([- ]*
>> 0>.1<.%.)&(-&S) E';'C
>>>> (R>:+/&.:*:)@:([- ]*
>>>> 0>.1<.(+/@:* -...@% ] +/@:* ]))&(-&S) E'
>>>> +----+-----+----+----+
>>>> |rank|tm*sz|time|size|
>>>> +----+-----+----+----+
>>>> | 1 |2.63 |1.18|2.23|
>>>> +----+-----+----+----+
>>>> | 0 |1.00 |1.00|1.00|
>>>> +----+-----+----+----+
>>>>
>>>> -:/".&>'C (R>:+/&.:*:)@:([- ]*
>> 0>.1<.%.)&(-&S) E';'C
>>>> (R>:+/&.:*:)@:([- ]*
>>>> 0>.1<.(+/@:* -...@% ] +/@:* ]))&(-&S) E'
>>>> 1
>>>>
>>>>
>>>> R.E. Boss
>>>>
>>>>
>>>>> -----Oorspronkelijk bericht-----
>>>>> Van: [email protected] [mailto:programming-
>>>>> [email protected]] Namens Marshall Lochbaum
>>>>> Verzonden: donderdag 7 oktober 2010 22:44
>>>>> Aan: 'Programming forum'
>>>>> Onderwerp: Re: [Jprogramming] Puzzle: line-circle/sphere
>> intersection>>>
>>>>> Alright. Take the projection of the center onto the line segment
>>>>> using %. .
>>>>> This is not allowed to be more than 0 or less than 1, so use
>>>>> (0>.1<.]) to take the endpoints if it is. Then find this
>> point on the
>>>>> line again and see if it is inside the circle.
>>>>> C (R>:+/&.:*:)@:([- ]* 0>.1<.%.)&(-&S) E
>>>>>
>>>>> Not tested yet, but I will update on that...
>>>>> Marshall
>>>>>
>>>>> -----Original Message-----
>>>>> From: [email protected]
>>>>> [mailto:[email protected]] On Behalf Of
>> Henry Rich
>>>>> Sent: Thursday, October 07, 2010 7:58 AM
>>>>> To: Programming forum
>>>>> Subject: Re: [Jprogramming] Puzzle: line-circle/sphere
>> intersection>>>
>>>>> I intend that a segment wholly inside the circle/sphere to
>> be counted
>>>>> as intersecting. In other words, the circle/sphere is
>> solid and we
>>>>> want to know if the segment touches any of it.
>>>>>
>>>>> Henry Rich
>>>>>
>>>>> On 10/7/2010 7:52 AM, R.E. Boss wrote:
>>>>>> That was bad reading from my side.
>>>>>> The formulas should be:
>>>>>>
>>>>>> ((S-C)<R)*.((E-C)<R) -...@+. ((S-C)>R)*.((E-C)>R)
>>>>>>
>>>>>> D=: +&.*:/"1 (S,:E)-"1 C
>>>>>>
>>>>>> (D<R) -...@+.&(*./) D>R
>>>>>>
>>>>>> Apart from that, the mathematics is not correct either. I
>> will give
>>>>> it
>>>>>> another thought.
>>>>>>
>>>>>>
>>>>>> R.E. Boss
>>>>>>
>>>>>>
>>>>>>> -----Oorspronkelijk bericht-----
>>>>>>> Van: [email protected] [mailto:programming-
>>>>>>> [email protected]] Namens Bo Jacoby
>>>>>>> Verzonden: donderdag 7 oktober 2010 13:02
>>>>>>> Aan: Programming forum
>>>>>>> Onderwerp: Re: [Jprogramming] Puzzle: line-circle/sphere
>>>>> intersection
>>>>>>>
>>>>>>> Hello Boss.
>>>>>>>
>>>>>>> I do not understand your notation. Henry wrote "startpoint
>> S and
>>>>>>> endpoint E". You wrote "endpoints A and B", but you use E
>> in your
>>>>> pseudo
>>>>> code.
>>>>>>> Please explain.
>>>>>>>
>>>>>>> It seems to me that you must compute the distances (s)
>> from C to S,
>>>>>>> (e) from C to E and (p) from C to the closest point P, (the
>>>>>>> perihelion), and you must determine (b) if P is between S
>> and E.
>>>>>>> Intersection occurs if (s<R and e>R) or (s>R and
>> e<R) or (b and p<R
>>>>> and
>>>>> (s>R or e>R)).
>>>>>>>
>>>>>>> --- Den tors 7/10/10 skrev R.E. Boss<[email protected]>:
>>>>>>>
>>>>>>> Fra: R.E. Boss<[email protected]>
>>>>>>> Emne: Re: [Jprogramming] Puzzle: line-circle/sphere intersection
>>>>>>> Til: "'Programming forum'"<[email protected]>
>>>>>>> Dato: torsdag 7. oktober 2010 09.49
>>>>>>>
>>>>>>> Since the word 'intersect' is used, I assume you mean the
>> boundary>>> of
>>>>>>> the circle/sphere (and not the area), so that a line
>> segment which
>>>>> is
>>>>>>> completely inside is not intersecting. One point at the boundary
>>>>>>> means intersecting, I assume.
>>>>>>> So the segment is not intersecting if and only if it is
>> completely>>>>> inside or completely outside the circle/sphere.
>>>>>>>
>>>>>>>
>>>>>>> If the line segment is given by the endpoints A and B, the
>> pseudo>>>>> code is rather straightforward:
>>>>>>>
>>>>>>> ((A-E)<R)*.((B-E)<R) -...@+. ((A-E)>R)*.((B-E)>R)
>>>>>>>
>>>>>>> so ones gets
>>>>>>>
>>>>>>> D=: +&.*:/"1 (A,:B)-"1 E
>>>>>>>
>>>>>>> (D<R) -...@+.&(*./) D>R
>>>>>>>
>>>>>>> (untested)
>>>>>>>
>>>>>>>
>>>>>>> R.E. Boss
>>>>>>>
>>>>>>>
>>>>>>>> -----Oorspronkelijk bericht-----
>>>>>>>> Van: [email protected]
>> [mailto:programming-
>>>>>>>> [email protected]] Namens Henry Rich
>>>>>>>> Verzonden: woensdag 6 oktober 2010 22:54
>>>>>>>> Aan: Programming forum
>>>>>>>> Onderwerp: [Jprogramming] Puzzle: line-circle/sphere
>> intersection>>>>>>
>>>>>>>> Given circle/sphere with center C and radius R, and a
>> line-segment
>>>>>>>> with startpoint S and endpoint E, write the J code to
>> tell whether
>>>>>>>> the line-segment intersects the circle/sphere.
>>>>>>>>
>>>>>>>> R is an atom, the rest are lists with 2 or 3 atoms.
>>>>>>>>
>>>>>>>> This problem arises in collision detection for games and
>>>>> simulators,
>>>>>>>> or if you are trying to see whether a path intersects a round
>>>>> obstacle.
>>>>>>>>
>>>>>>>> I found a solution whose brevity surprised me.
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
>
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
