컨볼루션, 고속 푸리에 변환 및 다항식 (2022)

 (alvarorevuelta.com)
1P by GN⁺ | ★ favorite | 댓글 1개
  • 큰 차수의 다항식을 고등학교식으로 전개하면 모든 항 쌍을 곱해야 하므로 O(n²) 비용이 빠르게 병목이 됨
  • 다항식의 계수 벡터 곱셈은 이산 신호의 컨볼루션과 같으며, [2, 3, 4]와 [5, 6, 7]의 결과는 [10, 27, 52, 45, 28]임
  • DFT는 이산 신호를 주파수 영역으로 옮기고, FFT는 같은 변환을 O(n log n)에 계산해 큰 입력에서 차이를 만듦
  • 시간 영역의 컨볼루션은 주파수 영역의 원소별 곱셈으로 바뀌므로, FFT로 변환해 곱한 뒤 IFFT로 되돌리면 다항식 곱셈을 더 빠르게 처리할 수 있음
  • 작은 차수에서는 FFT/IFFT 왕복 비용이 이득을 상쇄할 수 있지만, 차수가 커질수록 FFT 방식이 더 효율적임

다항식 곱셈이 느려지는 이유

  • 다항식 P(x)는 계수 a_k와 변수 x의 거듭제곱 항을 더한 형태로 표현됨
    • 예시 P(x)=5x²+2x+9는 차수가 2인 다항식임
    • 계수 벡터는 표기 방식에 따라 [5, 2, 9] 또는 [9, 2, 5]처럼 나타낼 수 있음
  • 덧셈과 뺄셈은 같은 차수의 항끼리 더하거나 빼면 되므로 비교적 단순함
    • Python에서는 zip(p, q)로 각 계수를 순회하며 a + b 또는 a - b를 계산할 수 있음
    • 차수가 다르면 zip_longest를 사용할 수 있음
  • 곱셈은 각 항을 서로 곱한 뒤 같은 차수 항을 다시 합쳐야 해 계산량이 커짐
    • (2x²+3x+4) × (5x²+6x+7)의 결과는 10x⁴+27x³+52x²+45x+28임
    • 이 방식의 복잡도는 O(n²) 이며, 차수가 커질수록 필요한 곱셈 수가 늘어남

계수 벡터와 컨볼루션

  • 이산 영역에서 두 신호 p와 q의 컨볼루션은 y[n]=Σ p[k]·q[n-k]로 정의됨
  • 계산은 q를 뒤집은 뒤 p 위에서 왼쪽에서 오른쪽으로 이동시키며, 겹치는 원소의 곱을 더하는 방식임
  • 예시 신호는 다음과 같음
    • p = [2, 3, 4]
    • q = [5, 6, 7]
  • q를 뒤집어 이동하면 각 출력 계수는 다음 순서로 만들어짐
    • 2×5 = 10
    • 2×6 + 3×5 = 27
    • 2×7 + 3×6 + 4×5 = 52
    • 3×7 + 4×6 = 45
    • 4×7 = 28
  • 컨볼루션 결과는 y = [10, 27, 52, 45, 28]임
    • 이는 다항식 곱셈으로 얻은 10x⁴+27x³+52x²+45x+28의 계수와 같음
    • 따라서 다항식 곱셈은 계수 벡터의 컨볼루션으로 볼 수 있음

푸리에 변환과 FFT

  • 푸리에 변환은 신호를 시간 영역에서 주파수 영역으로 변환함
    • 시간 관점에서는 특정 시점의 값으로 신호를 봄
    • 주파수 관점에서는 서로 다른 진동 주파수들의 합으로 신호를 해석함
  • 진동 주파수는 사인과 코사인으로 표현되며, 각각 계수와 위상을 가짐
  • 5Hz 순수 사인파에 FFT를 적용하면 주파수 영역에서 5Hz 위치에 델타처럼 나타남
    • 이는 시간 영역의 사인파가 5Hz 사인 하나로 표현될 수 있음을 보여줌
  • 관련 용어는 다음처럼 구분됨
    • Fourier Transform(FT): 연속 영역에서 정의된 푸리에 변환
    • Discrete Fourier Transform(DFT): 이산 신호에 대해 정의된 푸리에 변환
    • Fast Fourier Transform(FFT): DFT를 O(n²) 대신 O(n log n) 에 계산하는 알고리듬
  • DFT는 이산 시간 신호 x[n]을 주파수 영역의 X[k]로 바꿈
    • 각 X[k]는 입력 샘플을 특정 주파수를 나타내는 복소수와 곱해 더해 계산됨

주파수 영역에서 곱셈으로 바꾸기

  • DFT와 주파수 영역의 핵심 장점은 컨볼루션을 원소별 곱셈으로 바꿀 수 있다는 점임
    • 시간 영역에서 두 신호를 컨볼루션하는 것은 주파수 영역에서 두 신호를 곱하는 것과 같음
    • 곱셈은 컨볼루션보다 더 빠르게 계산할 수 있음
  • 다항식 곱셈을 빠르게 수행하는 절차는 다음과 같음
    • 다항식을 FFT로 주파수 영역으로 변환함: O(n log n)
    • 주파수 영역에서 원소별로 곱함: O(n)
    • 결과를 IFFT로 다시 시간 영역으로 변환함: O(n log n)
  • 전체적으로 FFT를 이용하면 다항식 곱셈을 O(n log n) 복잡도로 수행할 수 있음
  • 큰 다항식에서는 고등학교식 O(n²) 곱셈보다 빠름

Python 구현과 벤치마크

  • multiply_naive는 이중 반복문으로 모든 계수 쌍을 곱해 결과 위치 i + j에 더함
    • 결과 길이는 len(p) + len(q) - 1임
    • 복잡도는 O(n²) 임
  • multiply_fft는 FFT/IFFT 기반으로 계수 곱셈을 수행함
    • 결과 길이를 담을 수 있도록 len(p) + len(q) - 1 이상인 2의 거듭제곱 길이를 계산함
    • np.pad로 두 입력을 패딩함
    • np.fft.fft로 변환한 값을 원소별로 곱함
    • np.fft.ifft로 되돌린 뒤 실수부를 반올림해 정수 계수로 변환함
  • 예시 입력 p = [2, 3, 4], q = [5, 6, 7]에서는 두 방식 모두 [10, 27, 52, 45, 28]을 반환함
  • 벤치마크에서는 multiply_naive 대신 np.convolve를 사용하는 multiply_convolve와 FFT 방식을 비교함
    • multiply_naive는 Python 반복문이 느려, np.fft.fft를 쓰는 FFT 방식과 직접 비교하기 어렵기 때문임
    • np.convolve는 같은 연산을 저수준 C 코드로 수행함
  • 차수는 range(1, 30000, 1000) 범위로 늘리고, 각 차수에서 1부터 999999 사이의 임의 계수로 두 다항식을 생성함
    • 각 방식은 n_runs = 5로 평균 시간을 측정함
    • 낮은 차수에서는 FFT/IFFT 왕복 변환 비용 때문에 FFT 방식이 유리하지 않을 수 있음
    • 차수가 증가하면 FFT 방식이 훨씬 효율적인 결과를 보임

함께 보면 좋은 글 β

GN⁺   [-]
Hacker News 의견들

  • 이런 설명에서 늘 거슬리는 건 보통 수치 오차를 잊는다는 점임
    계수 곱셈을 그냥 “상수 시간”으로 추상화할 수는 없음. 그렇게 할 거면 애초에 곱셈 전체를 추상화해도 마찬가지임. 수치 정밀도를 고려하면 O(n (log n)^3)에 더 가까움 [1]
    [1]: http://numbers.computation.free.fr/Constants/Algorithms/fft....

    • 그 글에서 인용한 오차 한계는 지나치게 비관적임. Knuth 최신판에는 올바른 한계가 들어갔는데, 내가 알려줬기 때문임
    • OP 글에 나온 사원수 기반 연산을 활용하면 곱셈 오차를 줄이거나 아예 없앨 수도 있으면 좋겠음 [1],[2],[3]
      [1] One-Dimensional Quaternion Discrete Fourier Transform and an Approach to Its Fast Computation:
      https://www.mdpi.com/2079-9292/12/24/4974
      [2] Convolution Theorems for Quaternion Fourier Transform: Properties and Applications:
      https://onlinelibrary.wiley.com/doi/10.1155/2013/162769
      [3] On the Matrix Form of the Quaternion Fourier Transform and Quaternion Convolution:
      https://arxiv.org/abs/2307.01836
    • 계수가 정수라면 충분히 큰 법(modulus)을 쓰는 NTT로 정확한 결과를 얻을 수 있고, 특히 하드웨어에서는 곱셈 시간도 더 빨라질 수 있음
    • 그래서 컴퓨터 과학과 소프트웨어 엔지니어링을 구분하게 됨 :)

  • 이 방법으로 긴 숫자끼리 곱할 수 있음. 핵심은 다항식 곱셈이 올림(carry)을 하지 않은 일반 긴 숫자 곱셈과 같다는 점임
    예를 들어 1000자리 숫자가 있으면 각 자릿수를 1000개 원소를 가진 다항식의 계수로 삼음. 그런 다음 글에서 설명한 FFT 방식으로 이 다항식들을 곱할 수 있음. 결과를 다시 숫자로 바꾸려면 올림을 처리해야 함. 어떤 원소가 10보다 크면 초과분을 다음 자리로 넘기고, 계수들을 숫자로 바꾸면 됨
    기본 아이디어는 이렇고, 올림에 필요한 정밀도와 FFT 결과를 가장 가까운 정수로 반올림해도 맞는지 보장하는 부분에는 미묘한 점이 있음. 이런 방식이 이 분야 대표 라이브러리인 GMP에서 큰 수 곱셈을 하는 방식임

    • 말한 것처럼 10진수는 x=10인 다항식으로 표현할 수 있으니 이해됨. 예를 들어 983 = 9x^2 + 8x + 3, 즉 [9, 8, 3]임
      실제로 의미가 있으려면 숫자가 얼마나 커야 하는지, 그리고 어디에 쓰이는지가 궁금함

  • 아직 안 봤다면 이 영상을 보면 좋음
    https://youtu.be/h7apO7q16V0?si=bmgUEMTQSqU3flIv
    다항식 곱셈에서 FFT 알고리즘을 유도하는데 정말 훌륭함. 6개월쯤마다 다시 봄

  • FFT의 “합성곱은 점별 곱셈” 성질은 임의의 순환 곱셈군에서도 성립함. 더 대수적인 유도는 https://www.sciencedirect.com/science/article/pii/S002200007... 참고
    이를 “조화 FFT”라고 부르는 경우도 있고, 비조화 FFT도 있음: GF(2^n) 위의 [LCH14] “additive NTT”, 유한체의 단위원 X^2+Y^2=1 위의 [HLP24] circle FFT, 타원곡선 아이소제니열 위의 [BCKL21] ecfft
    [LCH14]: https://arxiv.org/abs/1404.3458
    [HLP24]: https://eprint.iacr.org/2024/278
    [BCKL21]: https://arxiv.org/pdf/2107.08473

  • 더 빠른 다항식 곱셈에 FFT를 쓰자고 처음 제안한 사람은 누구였을까?
    최근 궁금해서 찾아봤는데 인용 추적을 잘하진 못해도, David Eppstein의 1995년 논문 [0]까지는 거슬러 올라갔음. 여기서는 점진적 갱신 뒤 부분합 문제를 효율적으로 풀기 위해 이를 사용함. 분명 Knuth의 TAOCP에는 더 일찍 있었을 것 같음
    FFT 다항식 곱셈으로 반복 허용 정확한 부분합 문제를 준지수 시간에 풀 수 있다는 사실도 꽤 충격적이었음 [1]. 중요한 건 이 알고리즘이 O(N log N)인데, 여기서 N은 집합 크기가 아니라 최대 원소이므로 P ≠ NP에 대한 반례 같은 건 아님
    [0] https://escholarship.org/content/qt6sd695gn/qt6sd695gn.pdf
    [1] https://x.com/festivitymn/status/1788362552998580473?s=46&t=...

  • 모든 기계학습은 합성곱 방정식을 푸는 일이라고 생각함
    이 논문은 강화학습 맥락에서 다루지만 https://arxiv.org/abs/1712.06115, 대부분의 접근이 이 패러다임 안에 들어맞음

    • 기본적으로 커널 방법이라는 뜻 아닌가?

  • 방금 FFT를 사용해 시계열 부분열들의 큰 집합에 대한 내적을 계산하는 알고리즘(matrix profile)을 구현했음. 시계열 길이 n은 수억 단위까지 갈 수 있음
    FFT를 이용한 빠른 합성곱 계산으로 계산 시간이 O(n)에서 O(log n)으로 줄어들고, 이 규모에서는 속도 향상이 엄청남. GPU까지 쓰면 노트북에서 데이터 포인트 1천만 개를 0.1초에 처리하는 식으로 더 빨라짐

  • 이 연산의 핵심 “트릭”은 이 깨달음처럼 보임:

    다시 말해, 시간 영역에서 두 신호의 합성곱을 수행하는 것은 주파수 영역에서 두 신호를 곱하는 것과 같다.
    복잡한 아이디어를 훨씬 작은 단계로 나눠서, 수학이 약한 나도 어떻게든 이해할 수 있게 해준 좋은 글임. 하지만 중간 단계 하나를 놓친 걸까? 아니면 독자가 찾아보라는 연습문제로 남겨둔 걸까? 그 지점까지 이미 수학 능력을 한껏 끌어다 쓰고 있었는데, 약간 “그리고 나머지 빌어먹을 올빼미를 그리세요” 같은 느낌이었음. 나만 그런가? 글 자체는 정말 좋았음

    • 도움이 될지 모르겠지만: 학교에서 배우는 두 다항식의 곱셈은 실제로는 합성곱임
      “시간 영역에서 두 신호의 합성곱을 수행하는 것은 주파수 영역에서 두 신호를 곱하는 것과 같다”는 성질이 있고, FFT는 시간 영역에서 주파수 영역으로 변환할 수 있게 해줌. 그래서 다항식을 FFT로 주파수 영역으로 옮긴 뒤, 그 영역에서는 곱셈만 하면 됨. 합성곱보다 빠름. 빠진 단계가 이걸로 분명해지는지 궁금하고, 빠진 부분이 있으면 글을 업데이트할 수 있음

  • 그렇다면 정수 인수분해는 이산 디컨볼루션인가? FFT 표현, 즉 점별 곱셈의 역연산과 tableax, 즉 일반 긴 곱셈/올림 덧셈을 나란히 놓으면 대칭이 깨져서 빠른 알고리즘에 충분한 정보를 얻을 수 있을지 궁금함

  • 물론 순진한 다항식 곱셈은 다항식 차수에 대해 느림. 그런데 실제로 100차 다항식 두 개를 다뤄야 하는 때가 언제 있을까?
    이런 이유로 컴퓨터 대수 시스템에서는 이런 방법을 쓰지 않는다는 인상이 있음

    • 컴퓨터 대수 시스템, 예를 들어 Matlab의 chebfun은 임의의 함수를 100차 이상 다항식으로 바꿔서 근, 최적값 등을 더 쉽게 찾음
    • 오류 정정과 신호 처리에서는 매우 흔함
      https://www.youtube.com/watch?v=CcZf_7Fb4Us
      https://en.wikipedia.org/wiki/Reed%E2%80%93Solomon_error_cor...가 한 예임
    • 큰 파일의 CRC 체크섬 매개변수를 역공학하고 싶어서, 파일을 수백만 차 GF(2) 다항식으로 변환하고 최대공약수를 계산하는 프로그램[1]을 만들었음. FFT 기반 곱셈 없이는 합리적인 시간 안에 불가능함
      [1]: https://github.com/8051enthusiast/delsum
    • 이런 합성곱 관점과 FFT용 빠른 GPU 커널은 Mamba 이전의 일부 상태공간 모델에서 긴 시퀀스 모델링에 쓰였고, 여기서는 다항식이 입력 시퀀스임
      Hazy Research의 2020~2023년 블로그 글에 이 접근에 대한 정보가 많음
    • https://news.ycombinator.com/item?id=40306339 참고
      “(...) 물리학 연구에서 거의 1테라바이트 길이에 항이 1억 개가 넘는 표현식을 다룬 적이 있음”
답변달기