Алгоритм Зейделя. Haskell

[caxap]

Требуется реализовать алгоритм Зейделя для решения СЛАУ. Хочу реализовать на Хаскеле. Я Хаскель только учу и уверен, что я всё сделал не так, как сделали бы нормальные люди. Прошу поправить.

Про метод Зейделя (из Канатников, Крищенко "Линейная алгебра"):
Если имеется СЛАУ $Ax=b$, то мы сначала приводим её к нормальному виду $A^\top Ax=A^\top b$, т. к. чтобы алгоритм сходился, нужно, чтобы $A$ была положительно определённой и симметрической. $A^\top A$ такой является, и если исходная СЛАУ имеет решение, то и нормальная СЛАУ будет иметь такое же. Пускай СЛАУ $Ax=b$ уже нормальная. Тогда сам алгоритм:

$$x_i^{N+1}=\frac{1}{a_{ii}}\left(-\sum_{j=1}^{i-1} a_{ij} x_j^{N+1}-\sum_{j=i+1}^{n} a_{ij} x_j^N+b_i\right),\quad i=\overline{1,n}.$$

где $n$ -- число уравнений в СЛАУ. (В правой части есть $x_j^{N+1}$, но когда считается $x_i^{N+1}$, он уже подсчитан.)

(Код)

-- Тут можно не читать, я просто определяю некоторые матричные операции
type Vector = [Float]
type Matrix = [Vector]

infix 9 !
(!) :: Matrix -> (Int,Int) -> Float
m ! (i,j) = m !! (i-1) !! (j-1)

width :: Matrix -> Int
width m = length $ m !! 0

height :: Matrix -> Int
height = length

mul :: Matrix -> Matrix -> Matrix
mul a b | width a == height b =
    [[sum [a!(i,k) * b!(k,j) | k <- [1..n]] | j <- [1..n]] | i <- [1..m]]  where
        n = width a
        m = height b

transpose :: Matrix -> Matrix
transpose m = [[m!(i,j) | i <- [1..(height m)]] | j <- [1..(width m)]]

-- Вот сама функция решения СЛАУ
-- a', b' -- это матрицы A и b соответственно. n -- количество итераций (пока задаю сам).
solve :: Matrix -> Vector -> Int -> Vector
solve a' b' n | width a' == length b' = x n  where
    t = transpose a'    -- A^T
    a = t `mul` a'       -- A^T A
    b = (transpose $ mul t $ transpose [b']) !! 0      -- A^T b

    r = length b        -- n
    x 0 = replicate r 0
    x (n+1) = [(b!!(i-1) - sum [ a!(i,j) * (x (n+1))!!(j-1) | j <- [1..(i-1)]]
                - sum [ a!(i,j) * (x n)!!(j-1) | j <- [(i+1)..r]]) / a!(i,i)
              | i <- [1..r]]

Код работает:

-- {x+2y=4, 3x-4y=-2}, {x=1.2, y=1.4}
*Main> solve [[1,2],[3,-4]] [4,-2] 10
[1.1972656,1.3986328]

Но очень медленно (уже для $n=4$ и $N=8$ в ghci работает несколько секунд). Наверное я неправильно подошёл к реализации функции $x$. Она у меня сейчас возвращает список и, я боюсь, много раз вычисляется одно и то же. Была идея сделать бесконечный список решений $(x(0),x(1),\ldots)$, чтобы при вычислении $x(n+1)$ функция ссылаться к предыдущему элементу списка. Ничего не вышло

[.Serj.]

transpose есть в Data.List. Это раз.

Во-вторых, зря вы связались с length.

И, быть может, стоит вместо list comprehensions использовать монады или что-нибудь типа того?

//Тоже только учу Хаскелл.

P.S.

Цитата:

Она у меня сейчас возвращает список и, я боюсь, много раз вычисляется одно и то же.

http://www.haskell.org/haskellwiki/Memoization

[caxap]

.Serj.
Я знаю, что моя реализация матриц не очень хорошая. Но на фоне "слона" они на производительность не влияют. Матрицы маленькие и все матричные операции проводятся один раз: они нужны лишь для приведения СЛАУ к нормальному виду. На них можно вообще не смотреть (если СЛАУ изначально нормальная, то код с матрицами и половину solve можно выкинуть вообще). Основной код, повторяющийся много раз -- это

    x 0 = replicate r 0
    x (n+1) = [(b!!(i-1) - sum [ a!(i,j) * (x (n+1))!!(j-1) | j <- [1..(i-1)]]
                - sum [ a!(i,j) * (x n)!!(j-1) | j <- [(i+1)..r]]) / a!(i,i)
              | i <- [1..r]]

Про мемоизацию не понял. Ведь я даже не знаю, что он вычисляет много раз. Я просто "в лоб" записал формулу из книжки. Возможно не надо $x$ списком делать. А как тогда -- я не соображу...

[Xaositect]

caxap в сообщении #403002 писал(а):

А как тогда -- я не соображу...

    init = replicate r 0
    next xn = xn1
        where xn1 = [(b!!(i-1) - sum [ a!(i,j) * xn1!!(j-1) | j <- [1..(i-1)]]
                - sum [ a!(i,j) * xn!!(j-1) | j <- [(i+1)..r]]) / a!(i,i)
              | i <- [1..r]]

    iterations = init:(map next iterations)

-- Сб янв 22, 2011 14:34:27 --

Кроме того, я думаю, что использование более правильных структур данных для векторов и матриц улучшит производительность существенно. Посмотрите http://hackage.haskell.org/package/vector , а матрицу либо как Vector (Vector a), либо как обертку над вектором длины n*m (не помню, возможно там это делается автоматически с использованием type families)

[caxap]

Xaositect
Ого, летает. Спасибо! А почему старый вариант был такой тормознутый?

P. S. Что-то странно только. Маленькие системы хорошо решает, а дал из учебника -- решает неверно. (Старый код тоже не работает.)

*Main> solve [[7,-2,-3.2],[-22,5,-6],[7,8,9]] [0,1,-2] 100
[-6.270956e-2,-0.1374483,-5.127185e-2]
-- Тут нормально. Верный ответ: {{x -> -0.0627095, y -> -0.137448, z -> -0.0512719}}

*Main> solve [[0.1,-8,9.3,-8.2],[-3.3,2.4,-2.8,2.4],[-2.6,1.9,-2.2,1.9],[-1.3,9.3,-1.1,9.5]] [1.5,2,-1,0] 100
[-0.3682006,-0.11601713,0.13453485,7.84826e-2]
-- А тут не то. Верный ответ: {{x -> -224.851, y -> -12503.7, z -> 13.6381, t -> 12211.3}}

С чем это может быть связано?

[Xaositect]

caxap в сообщении #403039 писал(а):

Ого, летает. Спасибо! А почему старый вариант был такой тормознутый?

Думаю, выражения (x n) и (x (n+1)) вычислялись несколько раз.

caxap в сообщении #403039 писал(а):

С чем это может быть связано?

Не знаю. Я формулу не проверял, проверьте. Мб просто надо больше итераций, может матрица плохая.

[caxap]

Xaositect в сообщении #403045 писал(а):

Мб просто надо больше итераций, может матрица плохая.

Да нет, итерации явно сходятся к решению $\approx [-0.36821535,-0.11683835,0.13453573,7.928461e-2]$ (это для $N=500$). Самое интересное, что с многими другими СЛАУ работает. Формулу проверил много раз, она дословно переписана из учебника (формула из 1-го сообщения).
"Правильное" решение системы проверил в Mathematica. Причём как обычной, так и нормальной. Везде ответ тот же. То, как моя система находит нормальную СЛАУ, тоже проверил -- всё верно.

Мистика.

[Xaositect]

А еще у Вас в mul неправльно индексы расставлены, но вроде это для квадратных матриц не критично.

-- Сб янв 22, 2011 16:13:39 --

Хм
А в уравнение $A^{T}A x = A^{T}b$ это решение достаточно хорошо подходит.
Это матрица плохая. У нее число обусловленности огромное.

[caxap]

Xaositect в сообщении #403061 писал(а):

А еще у Вас в mul неправльно индексы расставлены

Ага, спасибо.

Xaositect в сообщении #403061 писал(а):

Это матрица плохая. У нее число обусловленности огромное.

Ааа. То есть авторы специально такую систему написали. А можно ли как-нибудь заставить находить алгоритм Зейделя настоящее решение? Или там надо будет очень-очень много итераций делать?

-- 22 янв 2011, 16:46 --

Пробовал менять $x_0$. Сходиться стало к другому, но всё равно не к тому, что надо.

[Xaositect]

Не, тут проблема не в том, что итераций надо много, а в том, что решение "скачет": небольшие изменения в матрице или векторе $b$ (а у нас неизбежно вылезут машинные погрешности) приводят к большим изменениям в решении.
Я из вычметодов уже мало что помню, спросите что ли в мат. разделе что делать если надо решить методом Зейделя плохо обусловленную систему.

[caxap]

Сделал функцию, чтобы с заданной точностью решала СЛАУ:

solve :: Matrix -> Vector -> Float -> Vector
solve a b eps = fst . head . filter (\(u,v) -> dist u v < eps) $ zip xs (tail xs) where
    xs = seidel a b  -- бывшая функция solve, возвращает беск. список приближений [Vector]

    norm :: Vector -> Float
    norm = sqrt . sum . map (^2)

    dist :: Vector -> Vector -> Float
    dist u v = norm $ zipWith (-) u v

Компилируется без ошибок. Но

*Main> solve [[1,1],[1,-1]] [5,2]

<interactive>:1:0:
    No instance for (Show (Float -> Vector))
      arising from a use of `print' at <interactive>:1:0-25
    Possible fix:
      add an instance declaration for (Show (Float -> Vector))
    In a stmt of an interactive GHCi command: print it

Никак не найду ошибку. Вроде бы с типами всё нормально. Я специально везде наставил определений через "::".

[Xaositect]

Так а где третий аргумент-то?
Компилятор, кстати, пишет, что не может вывести на экран функцию.

[Portnov]

Для хранения матриц лучше, наверное, использовать какой-нибудь вариант массивов вместо списков — см. модуль Data.Array и подобные, для большой производительности есть IOArray — полный аналог C'шных массивов. Также, надо хотя бы иметь ввиду, что для работы с матрицами (и, в частности, для решения СЛАУ) есть готовый пакет hmatrix (http://hackage.haskell.org/package/hmatrix), это привязки к LAPACK и GSL, так что производительность соответствует этим библиотекам.

[caxap]

Xaositect в сообщении #403145 писал(а):

Так а где третий аргумент-то?

Доучился...

Portnov
Учту. Хотя матрицы маленькие, они мало на производительность влияют.