473,320 Members | 1,883 Online
Bytes | Software Development & Data Engineering Community
Post Job

Home Posts Topics Members FAQ

Join Bytes to post your question to a community of 473,320 software developers and data experts.

Matlab: solving the linear system with QR/Householder

Hi all,
I'm trying to implement the QR method for solving the linear system Ax = b. The QR factorization is achieved using Householder method.

The main function is

Expand|Select|Wrap|Line Numbers
  1. function x = lin_solve(A,b)
  2. [R,v] = householder(A);
  3. y = Qt_times_b(v,b);
  4. x = R\y;
Here are the individual functions:

Expand|Select|Wrap|Line Numbers
  1. function [R,v] = householder(A)
  2. [m,n] = size(A);
  3. if m>=n,
  4.     NumberOfReflections = n;
  5. else
  6.     NumberOfReflections = m - 1;
  7. end
  8. R = A;
  9. v = cell(NumberOfReflections,1);
  10. for k = 1:NumberOfReflections,
  11.     x = R(k:m,k);
  12.     xnorm = norm(x);
  13.     if xnorm>0,
  14.         % Compute the normal vector of the reflector
  15.         v{k} = -x;
  16.         v{k}(1) = v{k}(1) - sign(x(1))*xnorm;
  17.         v{k} = (sqrt(2)/norm(v{k}))*v{k};
  18.         % Update R
  19.         for j = k:NumberOfReflections,
  20.                 R(k:m,j) = R(k:m,j) - (v{k}'*R(k:m,j))*v{k};
  21.         end
  22.     else
  23.         v{k} = zeros(m-k+1,1);
  24.     end
  25. end
  26.  
  27. function y = Qt_times_b(v,b)
  28. NumberOfReflections = length(v);
  29. y = b;
  30. p = NumberOfReflections;
  31. for k = 1:NumberOfReflections,
  32.    F = eye(p)-2*(v{k}*v{k}')/norm(v{k})^2;
  33.    p = p -1;
  34.    % Put F into "Q" to get Qk
  35.    y = Qk*y;
  36. end
I'm having trouble with the "Put F into Q to get Q_k" part.

I know Q is implicitly defined by the Householder reflectors stored in v, computed in the householder function, and y can be obtained by Q'b = Q_n * ... * Q_3 * Q_2 * Q_1 * b

Also, Qk =
[ I_(k-1) 0 ]
[ 0 F_k ]
matrix, where F_k = I - 2 (v * v') / norm(v)^2

But how do I put F into Q?

Thank you.

Regards,
Rayne
Sep 16 '07 #1
0 1989

Sign in to post your reply or Sign up for a free account.

Similar topics

2
by: Anatoly | last post by:
Hi All !!! Can anyone give me link to the C++ source or to the compiled program unit (compatible with Microsoft VC++), designed for the solving of the algebraic linear systems with rare matrices...
6
by: Lars Christiansen | last post by:
Hi! I am a master student in (geo)physics at the University of Copenhagen and part of a study group on C++ as a scientific programming language. I, and the other students in the group, have...
4
by: Thomas | last post by:
Hi, I've to solve a system of linear equations and inequations in C++. Are there (free) librarys for that topic? Regards, Thomas
0
by: m31hu1 | last post by:
Hello, i try to execute a .m file using PHP i wanted to count the number of jpg files in the directory but PHP was unable to detect it. the PHP script: <?php execute('start /b matlab...
6
by: omesh | last post by:
hi guys.. i got a problem.. could anyone help me write a program.. the program is as follows:- 2. Solving systems of Linear equations using Iteration: You are required to write a C++ program...
4
by: itcecsa | last post by:
Hi, I am implementing a small Python project, what I am going to do is to open Matlab and run some M-files, and get some output from Matlab command prompt. I have no idea how to open Matlab...
0
by: ursskmdali | last post by:
Hi, I have just created a sample in vc++ 6.0 to call MATLAB functions. Everything goes fine it compiles, builds fine with no errors and while executing it giving error like this "Missing ICU DATA...
5
by: adinda | last post by:
So what i need is this; (I'm very new at this,, programming in C I mean...) In matlab I had a while loop, and after each loop was done I added my resulting matrix to an object. Seeing the loop...
14
by: Luna Moon | last post by:
Dear all, Can C++/STL/Boost do the vectorized calculation as those in Matlab? For example, in the following code, what I really want to do is to send in a vector of u's. All other...
0
by: DolphinDB | last post by:
The formulas of 101 quantitative trading alphas used by WorldQuant were presented in the paper 101 Formulaic Alphas. However, some formulas are complex, leading to challenges in calculation. Take...
0
by: ryjfgjl | last post by:
ExcelToDatabase: batch import excel into database automatically...
0
by: Vimpel783 | last post by:
Hello! Guys, I found this code on the Internet, but I need to modify it a little. It works well, the problem is this: Data is sent from only one cell, in this case B5, but it is necessary that data...
0
by: jfyes | last post by:
As a hardware engineer, after seeing that CEIWEI recently released a new tool for Modbus RTU Over TCP/UDP filtering and monitoring, I actively went to its official website to take a look. It turned...
1
by: PapaRatzi | last post by:
Hello, I am teaching myself MS Access forms design and Visual Basic. I've created a table to capture a list of Top 30 singles and forms to capture new entries. The final step is a form (unbound)...
0
by: CloudSolutions | last post by:
Introduction: For many beginners and individual users, requiring a credit card and email registration may pose a barrier when starting to use cloud servers. However, some cloud server providers now...
0
by: Defcon1945 | last post by:
I'm trying to learn Python using Pycharm but import shutil doesn't work
0
by: Shællîpôpï 09 | last post by:
If u are using a keypad phone, how do u turn on JavaScript, to access features like WhatsApp, Facebook, Instagram....
0
by: af34tf | last post by:
Hi Guys, I have a domain whose name is BytesLimited.com, and I want to sell it. Does anyone know about platforms that allow me to list my domain in auction for free. Thank you

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.