# psykotic

By: a guest on Aug 9th, 2009  |  syntax: None  |  size: 17.54 KB  |  hits: 723  |  expires: Never
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
1. > import Data.List
3.
4. We start with some basic types and functions for vectors and matrices. Nothing there
5. is particularly interesting, so if you already know this material then you are encouraged
7.
8. ---
9.
10. Vectors can be scaled by a number and added in pairs when they have the same length.
11.
12. > vScale :: Num a => a -> Vector a -> Vector a
13. > vAdd   :: Num a => Vector a -> Vector a -> Vector a
14.
15. Vectors also have a dot product:
16.
17. > vDot :: Num a => Vector a -> Vector a -> a
18.
19. A vector will be represented as a list of numbers.
20.
21. > type Vector a = [a]
22.
23. Addition and scaling is then easily defined entrywise.
24.
25. > vScale a = map (a *)
26. > vAdd = zipWith (+)
27.
28. The dot product is computed by multiplying vectors entrywise and summing.
29.
30. > vDot v w = sum (zipWith (*) v w)
31.
32. Subtraction can be defined in terms of addition and scaling.
33.
34. > vSub v w = v `vAdd` ((-1) `vScale` w)
35.
36. A core operation in algorithms of linear algebra like Gaussian elimination is to
37. eliminate the first entry in a vector by subtracting an appropriate multiple of
38. another vector. This is only well-defined when the other vector's leading entry
39. is non-zero, or there will be a division by zero.
40.
41. > vEliminate (v:vs) (w:ws) = 0 : ws `vSub` ((w / v) `vScale` vs)
42.
43. ---
44.
45. A matrix is an mxn rectangular array of numbers. As with vectors, matrices can be
46. scaled by a number and added in pairs when equal in length.
47.
48. > mScale :: Num a => a -> Matrix a -> Matrix a
49. > mAdd :: Num a => Matrix a -> Matrix a -> Matrix a
50.
51. An mxn matrix can also be transposed to yield an nxm matrix.
52.
53. > mTranspose :: Num a => Matrix a -> Matrix a
54.
55. Unlike vectors, matrices can be multiplied, which is generally non-commutative.
56. An mxp matrix can be multiplied by a pxn matrix to yield a mxn matrix.
57.
58. > mMul :: Num a => Matrix a -> Matrix a -> Matrix a
59.
60. For every n, there is a special kind of nxn matrix, analogous to 1 for numbers,
61. which acts like a left and right identity for multiplication.
62.
63. > mIdentity :: Num a => Int -> Matrix a
64.
65. It is natural to think of a matrix in terms of its individual rows or columns.
66.
67. > mRows    :: Num a => Matrix a -> [Vector a]
68. > mColumns :: Num a => Matrix a -> [Vector a]
69.
70. A matrix will be represented as a list of row vectors.
71.
72. > type Matrix a = [Vector a]
73.
74. Its representation as a list of rows is then trivial.
75.
76. > mRows = id
77.
78. Scaling and adding is a simple matter of scaling and adding the rows.
79.
80. > mScale a = map (a `vScale`)
82.
83. To transpose a matrix, transpose the list of lists that is its representation.
84.
85. > mTranspose = transpose
86.
87. The columns of a matrix are the rows of the transposed matrix.
88.
89. > mColumns = mRows . mTranspose
90.
91. Multiplication is only a little trickier. The entries in the product consist of
92. every dot product of a row in the left operand with a column in the right operand.
93.
94. > mMul m1 m2 = [[row `vDot` column | row <- mRows m1] | column <- mColumns m2]
95.
96. Finally, a finite identity matrix is computed by taking a finite slice of an infinite
97. identity matrix.
98.
99. > mIdentity n = (map (take n) . take n) infiniteIdentity
100. >     where infiniteIdentity = iterate (0:) (1 : repeat 0)
101.
102. > mNumRows = length
103.
104. > mNumColumns []     = 0
105. > mNumColumns (x:xs) = length x
106.
107. TODO: Talk about the shape invariant. No empty rows, etc.
108.
109. ---
110.
111. The representation of a matrix as a list of rows automatically induces a data type of
112. matrix zippers that will turn out be very useful. In the sequel, for brevity 'zipper' will
113. without further qualification refer to matrix zippers rather than zippers generally.
114.
115. A matrix zipper is a  four-way block decomposition of a matrix that efficiently supports
117. the entries straddling the dividers. In diagrammatic form it looks like this:
118.
119.     / A | B \
120.     |---+---|
121.     \ C | D /
122.
123. The intersection between the dividers, marked by + in the diagram, is called the focus.
124.
125. Note that despite the symmetry of the diagram, the matrices A, B, C and D are not generally
126. all the same size. A and B have the same height, as do C and D. Likewise, A and C have the
127. same width, as do B and D.
128.
129. We will call A, B, C and D the upper-left, upper-right, lower-left and lower-right blocks
130. of the zipper, respectively.
131.
132. As the matrix representation is row-oriented, the induced zipper representation will
133. likewise be row-oriented, taking the form ((A, B), (C, D)).
134.
135. > type Zipper a = ((Matrix a, Matrix a), (Matrix a, Matrix a))
136.
137. One of the most important properties of zippers are their symmetries. The symmetry of a
138. plain old pair is swapping.
139.
140. > swap (a, b) = (b, a)
141.
142. A zipper has two levels of pairing and therefore two kinds of swapping. The first is vertical:
143.
144.     / A | B \     / C | D \
145.     |---+---| ==> |---+---|
146.     \ C | D /     \ A | B /
147.
148. > swapVertically = swap
149.
150. The second is horizontal:
151.
152.     / A | B \     / B | A \
153.     |---+---| ==> |---+---|
154.     \ C | D /     \ D | C /
155.
156. > swapHorizontally ((a, b), (c, d)) = ((b, a), (d, c))
157.
158. These swapping symmetries can be composed in two ways. The first is vertical then horizontal:
159.
160.     / A | B \     / C | D \     / D | C \
161.     |---+---| ==> |---+---| ==> |---+---|
162.     \ C | D /     \ A | B /     \ B | A /
163.
164. The second is the opposite, horizontal then vertical:
165.
166.     / A | B \     / B | A \     / D | C \
167.     |---+---| ==> |---+---| ==> |---+---|
168.     \ C | D /     \ D | C /     \ B | A /
169.
170. They give the same result!
171.
172. > swapDiagonally = swapVertically . swapHorizontally
173.
174. It's important to realize that though the diagram might suggest otherwise, 90-degree rotations are not
175. natural symmetries of zippers. In order to match up the invariants on width and height, we would need
176. to transpose the individual blocks:
177.
178.     / A | B \     / C^T | A^T \
179.     |---+---| ==> |---+-------|
180.     \ C | D /     \ D^T | B^T /
181.
182. This is a non-incremental (requiring O(n^2) time) transformation of the zipper's structure,
183. and therefore not considered a proper symmetry. In a sense this is the underlying reason for why a zipper
184. should not be considered exactly synonymous with a four-way block matrix.
185.
186. Thus a zipper has a total of four symmetries. The first is the identity, the trivial symmetry, which
187. does nothing. The next two are horizontal and vertical reflection. The final one is a rotational
188. symmetry that is associated with the diagonal and can be generated by composing horizontal and vertical
189. reflection in either order.
190.
191. Mathematicians know these symmetries as Klein's four-group. The underlying reason why the two orders
192. of composition for swapping are identical is that a 180-degree and a -180-degree rotation are the same.
193.
194. It was mentioned earlier that one of the purposes of zippers is to allow incremental motion of the dividers
195. and hence the focus. We will decompose motion into single steps along the the horizontal and vertical.
196.
197. Take the case of moving the vertical right by one step. With our representation of zippers this will entail
198. stripping off the column of B and adding it to A, and similarly for C and D. This is fundamentally similar to
199. moving the focus in a one-dimensional list zipper (A, B). There we would uncons a single entry from A and cons
200. that single entry onto B, assuming A is non-empty. Here our notion of unconsing and consing will have to be
201. more general because we are not dealing with single elements but with whole columns whose entries are scattered
202. throughout multiple rows.
203.
204. For this task a useful function will be unconsing a list into its head and tail and returning them as a pair:
205.
206. > uncons []     = Nothing
207. > uncons (x:xs) = Just (x, xs)
208.
209. We can extend this function over a matrix by unconsing all leading row elements and returning the remainder
210. rows as a matrix. When we uncons we have to ensure we normalize blocks so we don't leave behind any empty rows,
211. which would violate the invariant for matrix representation.
212.
213. > normalizeBlock ([]:m) = []
214. > normalizeBlock m      = m
215.
216. > unconsMatrix [] = Nothing
217. > unconsMatrix m  = do (v, m') <- mapAndUnzipM uncons m
218. >                      return (v, normalizeBlock m')
219.
220. The opposite of unconsMatrix is of course consMatrix, which conses a column vector onto a matrix. We add
221. an infinite supply of virtual empty rows onto the matrix to make sure no column entries are lost. Keeping
222. in mind what was said about block normalization for unconsing, this can be considered block denormalization.
223.
224. > denormalizeBlock = (++ repeat [])
225.
226. > consMatrix v m = zipWith (:) v (denormalizeBlock m)
227.
228. To step (A, B) right, we first check to make sure B is not empty. If it is then we do nothing. Otherwise we
229. uncons from B and cons onto A.
230.
231. > goRight' (a, b) = case unconsMatrix b of
232. >                   Just (v, b') -> (consMatrix v a, b')
233. >                   Nothing      -> (a, b)
234.
235. This is now easily extended to zippers pairwise.
236.
237. > goRight ((a, b), (c, d)) = (goRight' (a, b), goRight' (c, d))
238.
239. What about goLeft, goDown and goUp? Do we have to take pains and code them separately? This is where the
240. symmetry group of the zipper comes into play. Let's represent a right move as a diagram.
241.
242.     / A a | b B \     / A a b | B \
243.     |-----+-----| ==> |-------+---|
244.     \ C c | d D /     \ C c d | D /
245.
246. What happens if we first swap horizontally and then go right?
247.
248.     / A a | b B \     / B b | a A \     / B b a | A \
249.     |-----+-----| ==> |-----+-----| ==> |-------+---|
250.     \ C c | d D /     \ D d | c C /     \ D d c | C /
251.
252. This is a horizontal mirror image of a left move! Swapping horizontally again undoes the mirroring.
253.
254.     / A a | b B \     / B b | a A \     / B b a | A \     / A | a b B \
255.     |-----+-----| ==> |-----+-----| ==> |-------+---| ==> |---+-------|
256.     \ C c | d D /     \ D d | c C /     \ D d c | C /     \ C | c d D /
257.
258. This trick of applying a symmetry, then a function, then an inverse symmetry is called conjugation.
259. In the most general case the symmetry and its inverse will be distinct. Fortunately, all our symmetries
260. will be involutionary, meaning they are their own inverses.
261.
262. > conjugate g f = g . f . g
263.
264. > flipHorizontally = conjugate swapHorizontally
265. > flipVertically   = conjugate swapVertically
266. > flipDiagonally   = conjugate swapDiagonally
267.
269.
270. > goLeft = flipHorizontally goRight
271.
272. All the different motions of a zipper are generated by these single-step motions in cardinal directions.
273. That conjugating motion by symmetries gives another kind of motion (as opposed to some function that is a
274. non-motion) can be expressed mathematically by saying that the subgroup of motions is isotropic (or "normal")
275. within the zipper's group of automorphisms inside of which the zipper's symmetries reside.
276.
277. That takes care of horizontal motions. What about the vertical motions? Can we get them from horizontal motions
278. by conjugation? The only thing we haven't tried yet is conjugation by vertical swapping. Unfortunately horizontal
279. motions and vertically swapping commute for much the same reason that horizontal and vertical motions commute.
280. Not too surprising when you think about it. The issue is that the subgroup of horizontal motions and the subgroup
281. of vertical motions are individually isotropic.
282.
283. We are in fact stuck. Earlier we talked about how 90-degree rotations are not real zipper symmetries. If they
284. were, we could easily express goDown as the conjugation of goRight by a clockwise 90-degree rotation. Here the
285. gods of mathematics led us down the right path by excluding these as symmetries. If we had implemented goDown
286. in terms of 90-degree conjugation of goRight then it would have taken O(n^2) rather than O(n) time because of
287. the O(n^2) cost of matrix transposition.
288.
289. Fortunately there is a simple direct and efficient O(n) implementation of goDown. Let's look at a diagram.
290.
291.     / A | B \     / A | B \
292.     |---+---|     | c | b |
293.     | c | d | ==> |---+---|
294.     \ C | D /     \ C | D /
295.
296. When we implemented goRight, we only had to uncons the first column, which took a tiny bit of work. Here we
297. only have to uncons the first row from C and cons it onto A. It's exactly like shifting the focus in a
298. one-dimensional list-based zipper.
299.
300. > goDown' (a, b) = case uncons b of
301. >                   Just (v, b') -> (v : a, b')
302. >                   Nothing      -> (a, b)
303.
304. As with goRight', it is easy to extend this to full zippers.
305.
306. > goDown ((a, b), (c, d)) = ((a', b'), (c', d'))
307. >     where (a', c') = goDown' (a, c)
308. >           (b', d') = goDown' (b, d)
309.
311.
312. > goUp = flipVertically goDown
313.
314. That takes care of most of the groundwork we needed to lay down. Only a few more things remain.
315.
316. We've already seen how important the vertical and horizontal symmetries are in understanding zippers.
317. Each symmetry axis is associated with an interesting kind of structure, something close to a bialgebra.
318.
319. The usual kind of multiplication is something that takes two elements and produces a new element:
320.
321.     (a, a) -> a
322.
323. What happens if you turn the arrow around?
324.
325.     a -> (a, a)
326.
327. You get a comultiplication! Here we are of course dealing with zippers, so it takes this form:
328.
329.     Zipper a -> (Zipper a, Zipper a)
330.
331. What's a way to produce two zippers from a zipper? We cut it two! There are two axes and therefore
332. two comultiplications. For the horizontal axis, this is what it looks like:
333.
334.                     / A | B \
335.     / A | B \       \---+---/
336.     |---+---|  ==>
337.     \ C | D /       /---+---\
338.                     \ C | D /
339.
340. It looks simple and is even simpler to implement.
341.
342. > splitHorizontally ((a, b), (c, d)) = (((a, b), ([], [])), (([], []), (c, d)))
343.
344. Vertical splitting is no more difficult.
345.
346.     / A | B \       / A |\  /| B \
347.     |---+---|  ==>  |---+|  |+---|
348.     \ C | D /       \ C |/  \| D /
349.
350. > splitVertically ((a, b), (c, d)) = (((a, []), (c, [])), (([], b), ([], d)))
351.
352. What's the natural kind of multiplications associated with these comultiplications? It
353. seems pretty obvious. If you read the diagrams in reverse and it looks like concatenation.
354. The only hitch is that this assumes the focuses (the + marks in the diagrams) match up.
355.
356. If we want to define multiplication in this way for arbitrary zippers (at least of compatible
357. dimensions) then we will need to somehow normalize the focus before trying to concatenate.
358. It's clear there are only four possible distinguished focus positions for a zipper: the four corners.
359.
360. Moving the focus to a given corner is most easily accomplished by moving diagonally until an
361. empty block in that direction is reached. To go to the upper-right corner, for example,
362. we would iterate goRight . goUp until we stop making progress. For this a simple utility
363. function will be useful.
364.
365. > iterateUntil p f x = if p x then x else (iterateUntil p f) (f x)
366.
367. It turns out we need only define one of the four functions. We pick goUpperLeft.
368.
369. > goUpperLeft  = iterateUntil isCorner (goUp . goLeft)
370. >     where isCorner ((a, b), (c, d)) = all null [a, b, c]
371.
372. Earlier we discovered to our delight that we could get goLeft for free from goRight by
373. symmetric conjugations. We were not quite so delighted to discover that we couldn't get the
374. vertical motions this way, so we were forced to manually define at least one of the vertical motions.
375. However, in the present case the direction of motion for goUpperLeft is diagonal, and a
376. diagonal vector is non-invariant under both vertical and horizontal reflections. These
377. reflections generate all four 90-degree rotations of the diagonal motion!
378.
379. > goUpperRight = conjugate swapHorizontally goUpperLeft
380. > goLowerRight = conjugate swapVertically   goUpperRight
381. > goLowerLeft  = conjugate swapHorizontally goLowerRight
382.
383. With these functions, we can now easily implement conversion between matrices and zippers.
384. Conversion from a matrix to a zipper happens by putting the matrix in D's place and filling
385. A, B and C with blanks.
386.
387.      /+---\
388.      || M |
389.      \|   /
390.
391. > toZipper m = (([], []), ([], m))
392.
393. Converting a zipper to matrix conceptually erases the dividers between A, B, C and D.
394.
395.     / A B \
396.     \ C D /
397.
398. Equivalently, it moves the focus to the upper-left corner and returns the lower-right block.
399.
400.     /+-----\
401.     || A B |
402.     \| C D /
403.
404. > fromZipper z = let (([], []), ([], m)) = goUpperLeft z in m
405.
406. That was a small but crucial detour. Now we can get back to splitting and concatenating.
407.
408. To horizontally concatenate two zippers of the same height, we will have to move the
409. focus of the left zipper one of its two right corners, and conversely for the right zipper.
410. That leaves a choice of a top and bottom corner. Though a consistent choice of either
411. would yield a well-defined concatenation, the upper-left corner of a zipper is in a sense
412. the most canonical corner. That is to say, we are really working with pointed zippers, those
413. that have a distinguished origin. Thus we will always pick top over bottom and left over right
414. whenever there is a choice like that to be made.
415.
416. With those preliminaries out of the way, concatenation is surprisingly easy to implement.
417.
418. > concatHorizontally zLeft zRight = (([], []), (c, d))
419. >     where (([], []), (c, [])) = goUpperRight zLeft
420. >           (([], []), ([], d)) = goUpperLeft  zRight
421.
422. > concatVertically zTop zBottom = (([], b), ([], d))
423. >     where (([], b), ([], [])) = goLowerLeft zTop
424. >           (([], []), ([], d)) = goUpperLeft zBottom
425.
426. TODO: More about bialgebra properties that intertwine concatenation and splitting.
427.
428. ---
429.
430. > iterateTimes 0 f = id
431. > iterateTimes n f = (iterateTimes (n-1) f) . f
432.
433. TODO: Gauss, Laplace, Gauss-Seidel, etc.